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SYNOPSIS 

RAM SAJIWAN GUPTA 
Ph.D. 

Department of Mechanical Engineering 
Indian Institute of Technology Kanpur 
August , 1 977 

AUTOMATED OPTIMUM DESIGN OF AXIAL FLOW GAS TURBINE STAGE 

A computational capability for the automated optimum 
design of axial flow gas turbine stage to satisfy aerodynamic 
and mechanical constraints is developed in the present worh. 

More specifically an axial flow gas turbine stage is designed 
for maximum efficiency and/or minimum weight to satisfy 
strength, natural frequency and deflection requirements along' 
with the behavior constraints due to aerodynamic considerations. 

A survey of the available literature shows that 
mostly gas turbines are designed by conventional methods either 
from aerodynamic or mechanical considerations and the use of 
computers in this area is confined to the study of flow 
behavior ^or performance analysis. The study of vibration and 
strength has been carried out for idealized individual blades 
without considering their interaction with actual turbine 
design. The potentialities of optimization techniques have 
not yet been fully exploited by the turbomachine designers . 

In the present work an attempt has been made to develop a 
unified procedure for the design of axial flow gas turbines 



XX 


and to replace the conventional turbine design processes by more 
sophisticated design methods developed recently. 

In the present work the mean diameter of turbine , rotor 
blade chord to mean diameter ratio, nozzle blade chord to mean 
diameter ratio, nozzle blade spacing to mean diameter ratio, 
rotor blade spacing to mean diameter ratio, gas angle relative 
to the rotor blade at inlet, gas angle relative to the rotor 
blade at outlet and axial velocity of flow across the stage are 
taken as the design variables. 

A linear combination of losses (one minus efficiency) 
and weight (mass) of the gas turbine stage is taken as the 
objective for minimization. The major losses like profile loss, 
tip clearance loss and secondary loss have been accounted for 
in calculating the efficiency. The three dimensional effects 
have been included by assuming a free vortex flow . The effect 
of Reynolds number on the stage efficiency has also been 
considered. The various air properties used in the computations 
are obtained by fitting polynomial equations in the data using 
least squares method . The weight (mass ) of the stage is taken 
as the sum of the weights (masses) of the disc, rotor blades 
gnd nozzle blades . An iterative method has been developed for 
,^e computation of the efficiency of the stage and the method 
has been found to converge within three to four iterations. 
Bounds are placed on the flow coefficient, stage temperature 
drop coefficient, degree of reaction at the mean radius and at 
the root, and Mach number by considering the aerodynamic 
behavior of the gas flow in the turbine stage. 
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The vibration analysis of beams ha.ving rectangular 
cross sections has been made using finite element technique. 

A new beam element has been developed by considering the 
bending deflection, bending slope, shear deflection and shear 
slope in two mutually perpendicular planes as nodal degrees of 
freedom. V/ith these degrees of freedom the coupling effect due 
to pretwist can be accoimted for accurately. The effects of 
taper, twist, rotation, shear deformation and rotary inertia 
are considered in deriving the stiffness and mass matrices of 
the element. The various special cases from this general' 
element are also derived. The first four natural frequencies 
and mode shapes have been computed for several cases and the 
effects of breadth and depth taper ratios, twist angle, shear 
deformation, offset and rotation on the natural frequencies of 
vibration of cantilever beams are studied. The present 
numerical results are compared with those available in the 
literature wherever possible and have been found to be in good 
agreement even when the beam is idealized with only four finite 
elements . 

The tapered and twisted blade of airfoil cross section 
has been converted into an equivalent tapered and twisted 
cantilever beam of rectangular section for the purpose of 
vibration and stress analysis. The forces due to pressure 
and gas bending are assumed to act at the nodal points of the 
finite elements. The tip deflection and the stresses due to 
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the nodal forces are computed before calculating the natural 
frequencies of the idealized blade. The centrifugal stress, 
which can be calculated from the known dimensions of blade, is 
also included in computing the total stress induced in the 
blade. The Cholesky decomposition of symmetric banded matrices 
has been applied in solving the equilibrium equations and in 
obtaining a partial solution to the eigen value problem by 
Rayleigh -Eitz sub-space iteration algorithm. This algorithm 
solves the eigen value problem directly without a transformation 
to the standard form. 

The constrained optimum design problem is cast as a 
nonlinear mathematical programming problem. The interior 
penalty function method, with a variable metric unconstrained 
minimization technique is used to solve the optimum design 
problem. The computation of the gradient and slope of the 
function has been carried out by backward finite difference 
method . 

Five problems have been solved to demonstrate the 
effectiveness and generality of the optimization procedure 
developed. Out of these one problem has been solved by consi- 
dering only side constraints in order to find the effect of the 
behavior constraints on the optimum point. The results of 
this problem show that a lesser value of objective can be 
obtained by neglecting the behavior constraints; however, some 
of the behavior (response) quantities may take abnormal values 
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at the optinun point as no bounds are placed on them. It has 
also been found that the effects of rotation, shear deformation 
and rotary inertia do not change the optimum results appreci- 
ably. In most of the cases, the stress at the root and the 
degree of reaction at the root of the rotor blade have reached 
their limiting values at the optimum point. Two different 
starting points taken in the case of maximization of efficiency 
gave approximately the same results. A sensitivity anaJ.ysis 
has been conducted to find the effect of each of the design 
variables on the objective as well as on the response 
quantities. It has been found that efficiency is more sensitive 
to the chord and the spacing of the rotor blades while weight 
(mass) is most sensitive to the mean diameter of the gas 
turbine rotor. 

The contributions of the present work can be summarized 
as follows: (i) Development of a beam finite element for the 
vibration analysis of rotating, doubly tapered and pre-twisted 
Timoshenko beams, (ii) Application of the finite element 
developed for the deflection, stress and vibration analysis of 
rotor blades having airfoil section, (iii) Computerization of 
air properties for application in gas flow design problems, 

(iv) Development of a computer programme for the automated 
optimum design of gas turbine rotor stage with aerodynamic and • 
mecharacal constraints, (v) Study of sensitivity of losses, 
weight and other response quantities of the turbine stage with 
respect to the design variables. 



CHAPTER 1 


INTRODUCTION 


The application of classical aerodynamics in the design 
of turbomachinery, made the gas turbine to become a reality 
during the second world war. By replacing the conventional 
power units, gas turbines are now a days being being used 
successfully in space, aviation, marine, nuclear, rail -road, 
vehicular, cryogenic and petro -chemical applications. The 
rapid application of gas turbine plants in all phases of power 
generation is credited not only to the many advantages inherent 
in these plants, but also to the current technological progress. 
In gas turbine designs, the attainment of high temperature and 
high rotative speed alongwith the associated high efficiency 
has been possible due to the metallurgical advances made and the 
background of fluid mechanics, thermodynamics and gas dynamics 
accumulated . 

The analytical procedures used in designing turbomachines 
have now undergone a major revision. The empirical relations 
used in the design process are being replaced by more sophisti- 
cated methods of analysis. These methods provide detailed 
information regarding static and dynamic loads as well as 

machine response in the form of stresses and deflections. | 

i 

I 

In recent years, turbines of considerable size and ! 

i 

capacity are being produced wilh many stages. Along with the | 
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demand for larger gas turbines, factors like efficiency and 
component weight have become major design considerations. The 
reduction in weight , while maintaining safe stress levels , has 
required a more detailed knowledge of load and stress distribu- 
tion. The increased component flexibility resulting from weight 
reduction will lead to larger deflections that must be controlled 
during steady state and transient operating conditions. 

In the present work an attempt has been made to maximize 
the efficiency and/or minimize the weight of an axial flow gas 
turbine stage by considering deflection, stress and vibration 
aspects along with aerodynamic requirements. 

1 . 1 Review of literature . 

Since the present work deals with the automated optimiAin 
design of gas turbine rotor stage with mechanical design consi- 
derations, the available literature is reviewed under the head- I 
ings of aerodynamic design of gas turbines, mechanical design of | 
gas turbine blades and optimum design of turbomachines. | 

I 

1.1.1 Aerodynamic Design of Gas Turbines I 

The major developments in the design of gas turbines | 

began during the second world war. In his first authoritative | 

1 

book on the subject of steam and gas turbines, Stodola has 

i 

given the details of earliest practical gas turbine developed | 

2 ' 

by Holzwarth in 1908. In 1948, Reeman published the details | 
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of a simple jet propulsion turbine for aircraft application. 
Emmert"' gave a procedure for designing gas turbines for power 
plants using two-dimensional theory by considering mechanical 
aspects, Georgian^ published the design data of an experimental 
gas turbine used at Bengal Engineering College, Shibpur, India. 

5 

Carter discussed various approaches for incorporating 

three dimensional design modifications. Johnston and Knight^ 

reported the comparative results of two and three dimensional 

designs having free vortex blading for a single stage gas 

turbine of radius ratio 1.37. Wu and Wolfenstein applied 

streamline curvature method to satisfy the radial equilibrium 

conditions in axial flow compressor and turbine designs . 

8 

Hawthorne and Ringrose developed the actuator disc theory of 

compressible flow in free vortex turbomachinery to take into 

account the three dimensional effects in turbomachine design. 

g 

In his book, Horlock has reported the plots of Carmichael for 

determining the axial velocity distribution in the constant angle ‘ 

10 

method of three dimensional design. Smith and Johnston | 

designed, fabricated and conducted the performance tests of an i 

i 

experimental single stage turbine and also studied the appli- | 

1 1 

cability of Wu's method of design. Shaw presented a method i 

of designing blades in which the design of blade roots for i 

1 2 

heavy duty operation was discussed. Davis and Millar | 

I 

compared the matrix and stream line curvature methods for axial I 
flow turbomachines. A number of works have been reported at the ' 
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symposium entitled "Technical Advances in Gas Turbine Design" 

1 3 

conducted by the Institution of Mechanical Engineers . The 

NASA special publications^"^’ give comprehensive references in 

the area of gas turbine analysis, design and application. 

The application of high speed digital computers is also 

being made in solving design analysis and flow problems of gas 

turbines. Carter, Platt and Denherr developed a computer 

program, analysed the geometry and design point performance of 

axial flow turbines by making use of stream -filament procedure 

and developed a correlation for total pressure loss coefficient. 
17 

Glassman also reported the development of a computer programme 
for predicting the design analysis of axial flow turbines. 
Wasserbauer and Glassman developed a programme for predicting 
the off -design performance of radial inflow turbines. 

1.1.2 Mechanical Design of Gas Turbine Blades 

The blade design is not complete without the consider- 

1 9 

ation of stress, deflection and vibrational aspects. Hodge 

considered the stress analysis of gas turbine blades. 

20 

Pollmann included temperature effects in analysing hollow gas 

21 

turbine blades . Manson studied stresses in the blades and 
discs of gas turbines. Most of these investigators used conven- 
tional methods of strength of materials in their analyses. 

Further information regarding stresses in blades has been 

22 23 

compiled by Cox and Sawyer in their handbooks on gas 
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turbines. Very little information is available on the magnitude 

of defu-ectipn in gas turbine blading. 

Considerable amount of work has been published on the 

vibration analysis of single turbomachinery blade. The important 

references in this direction are shown inl'fgurel.1 under the 

categories of tapered beam, rotating beam, pre-twisted beam, 

asymmetric beam, and uniform beam (by finite element method). 

The effect of taper on vibration has been considered by 

Thomson ^ using matrix methods. Martin ^ analysed the tapered 

beam vibration problem by using perturbation method. Carnegie 
26 

and Thomos obtained the natural frequencies of long tapered 

27 

blades using finite difference technique. Rao and Carnegie 

determined the frequencies of lateral vibration of tapered 

cantilever beams by the use of Ritz-G-alerkin process. Mabie and 
28 

Rogers obtained the frequencies of vibration of doubly tapered 
cantilever beams with end mass using Bessel functions and 
presented the results in dimensionless form for various taper 
ratios. 

The effect of rotation on the vibration of beam has been 

29 

considered by many investigators. Plunkett studied free and 

forced vibration of rotating blades using matrix methods. Bo 
30 

and Renbarger found the bending vibrations of a rotating beam 

by applying Rayleigh -Ritz method. Boyce, DiPrima and 
31 

Handelman reported about the vibration of rotating beams of 
constant cross section using Rayleigh-Ritz and Southwell 
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32 

procedure. Tntema presented a simplified procedure and 

charts for rapid estimation of bending frequencies of rotating 

beams using Rayleigh-Ritz method. Carnegie, Stirling and 

Rleming^^ studied the vibration characteristics of turbine 

blading using finite difference technique. 

The effect of pre-twist has been studied by many 

34 

research workers. Mendelson and Gendler used station func- 
tions in analysis and performed experiments to determine the 

35 

effect of twist on the vibration of cantilevers. Rosard 
found the natural frequencies of twisted cantilevers by 
Myklestad method. DiPrima and Handelman'^ , Carnegie^' and 

'ZO 

Dawson^ applied the Rayleigh-Ritz method for the vibration 

39 

analysis of twisted beams. Rao analysed the flexural 
vibrations of pretwisted beams of rectangular cross section by 
Galerkin procedure. 

The vibration of beams having asymmetrical section 

40 

has also been considered by many authors. Garland obtained 

the normal modes and frequencies of beams having noncollinear 

' 41 

elastic and mass axes using Rayleigh-Ritz method. Targoff 

presented the matrices associated with the bending and coupled 

bending -tors ion vibration of beams using Holzer-Myklestad 

42 

procedure. Mendelson and Gendler determined the coupled 

bending-torsion vibrations of cantilever beams by means of 

43 

station functions, Carnegie and Dawson obtained the 
vibration characteristics of straight blades of asymmetrical 
airfoil cross section by applying transformation method. 
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The combined effects of taper, pre-twist, rotation and 

4-4 

assrmmelxy has also been reported in the literature. Targoff 

considered twist and rotation to study the bending vibrations 

45 

of beams using matrix methods. Rao and Carnegie presented a 

numerical procedure for the determination of frequencies and 

mode shapes of lateral vibration of blades allowing for the 

effects of pre-twist and rotation. Houbolt and Brooks solved 

the' differential equations of motion for combined flapwise 

bending, chordwise bending and torsion of twisted nonuniform 

rotor blades by using Rayleigh -Ritz method. Carnegie, Dawson 
47 

and Thomas applied finite difference technique to study the 
vibration characteristics of cantilever blading including the 
effects of asymmetry and pre -twist. Carnegie and Thomas^ 
obtained the coupled bending -bending frequencies of pre -twisted 

tapered blading using finite difference technique. Banerjee : 

49 ! 

and Rao applied Galerkin method to analyse the coupled I 

bending -tors ion vibrations of rotating blades and compared their t 

50 ^ 

results with experimental findings. Isakson and Eisley and i 

51 ‘‘ 

Montoya found the natural frequencies in coupled bending and 

torsion of twisted rotating and nonrotating blades. Krupka I 

52 ' 

and Baumanis included the effects of rotary inertia and ; 

i 

shear deflection in analysing the bending -bending mode of ; 

I 

vibration of a rotating tapered and twisted turbomachine blade | 
using Myklestad method. | 
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The finite element technique has also been applied by 

many ir vestigators , mostly for the vibration analysis of uniform 

beams. Ail these investigations differ one from the other in 

the nodal degrees of freedom taken for deriving the elemental 

53 

stiffness and mass matrices. McCalley derived the consistent 

mass and stiffness matrices by selecting total deflection and 

54- 

total slope as nodal coordinates. Archer analysed beams having 

55 

different boundary conditions. Kapur took bending deflection, 

shear deflection, bending slope and shear slope as nodal degrees 

of freedom and derived the elemental matrices of beams having 

56 

linearly varying inertia. Carnegie, Thomas and Documaki 

analysed uniform beams by taking few internal nodes in the 

57 

element. Nickel and Sector used total deflection, total 

slope and bending slopes of two ends and bending slope of the 

mid point of the element as degrees of freedom to derive the 

element stiffness and mass matrices of order seven. Thomas and 
5b 

Abbas analysed uniform Timoshenko beams by taking total 
deflection, total slope and derivatives of bending slope as 
nodal degrees of freedom. 

The study of banded group of blades has been reported 

59 

by some investigators. Campbell made an experimental investi- 
gation of low pressure blades grouped by lashing wires and 
stated an empirical relation between the natural frequencies of 
rotating and non-rotating blades. Prohl^^ analysed the problem 
of high pressure turbine blades by Holzer's technique and 
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obtained the first two natural frequencies. In a subsequent 

r -1 

paper I'rohl^ presented a case study for determining vibrational 

amplitude and stresses at resonance. Bhide^^ studied the free 

vibrations of packetted high pressure blades by assuming the 

blades to be straight, untwisted and uniform cantilevers and 

•65 

applied finite element method to get the solution. Baj^^; 
studied the free vibration characteristics of packetted turbine 
blades of high pressure section using finite element method by 
considering the blades to be untwisted tapered cantilevers . 

1 .1 .3 Optimum Design of Turbomachines 

Optimization techniques are versatile and can be 
applied to a large class of engineering problems. So far they 
have been mostly applied to structural engineering problems, 
and to a lesser extent, to other types of design problems. Box 
and Kapoor^"^ reported a capability for the minimum weight 
optimum design of planar truss -frame structures with inequality 
constraints on the maximum dynamic displacement , stress and 
natural frequencies using the method of feasible directions* 

Rao^^ presented the optimum design of aircraft wings with 
strength, stability, frequency and flutter constraints usir^ 
finite element idealization and interior penalty function 

approach. A review of the application of optimization techniques 

6 6 

in mechanical design has been published by Seireg . Raphael 
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and James published a programme through which a wing section 

68 

was op'.imized for weight. Reddy and Rao applied optimization 

techniques in the automated optimim design of machine tool 

structures with constraints on rigidity and chatter stability. 

6 9 70 

Hati and Rao ’ presented deterministic and probabilistic 

procedures for the determination of optimum machining conditions 

71 

for jobs involving single and multiple operations. Rao 

reported a capability for the minimum weight design of bridge 

girders for electric overhead travelling cranes. 

Very few papers are available regarding the application 

of optimization techniques in the design of turbomachines. 

72 

George reported about the optimization of a rocket engine 

73 

turbine using differential calculus approach. Balje applied 

optimization techniques in the design of axial flow turbines 

to obtain maximum stage efficiency. He solved the problem by 

considering six design variables; nozzle angle at the outlet 

of stator blade, blade angle at the outlet of rotor blade, blade 

height to rotor diameter ratio, blade chord to rotor diameter 

ratio of nozzle blades , blade chord tor rotor diameter ratio 

of rotor blade, and the degree of partial admission. Upper 

and lower bounds on the design variables were the constraints 

and Wood's technique, which is essentially a pattern search 

method, was employed in obtaining the solution of the problem. 

74 

Swift published a flow chart for optimizing pump-turbine 
designs using computers and reported the transient behavior of 
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75 

a pump -turbine . Kar and Reddy found the optimum shape of the 
impeller of a pump by using differential calculus approach. 

ry ^ 

Gupta' optimized the radial impeller of a pump by using 

graphical procedure. Saravanamuttoo and Macisaac ' published 

about the use of hybrid computer simulation of single -spool 

turbojet engine for optimizing the thrust response of gas 

78 ■ 

turbines. Yadav and Gupta optimized the channel circulation 

79 

loss of a centrifugal compressor. Blaho discussed the 
optimum design of axial flow fans from the view point of losses 

on 

using experimental curves . Paranjpe and Murthy discussed 
about the optimization and standardization of steam turbine 
blade profiles. In this work the authors discussed only a 
comparative procedure without giving any specific data of the 
problem. 

1.2 Objective and Scope of the Present Work : 

It can be seen from the available literature that the ; 

I 

potentialities of optimization techniques have not been exploited I 

I 

much by gas turbine designers. It is observed that the design- j 
ers considered either aerodynamic aspects only or mechanical i 

aspects only in designing gas turbines. This might be due to I 

i' 

the complex aerodynamic fluid flow behavior present in toirbines ! 
that gives rise to troublesome behavior constraints v/hich, in I 

the presence of constraints on the mechanical behavior, makes ' 

the problem more complicated . In the present work the automated - 
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optimum design of axial flow gas turbine stage is considered 
by including both aerodynamic and mechanical, considerations. 

The weight of the stage and the losses of the stage 
have been considered in formulating a compound objective 
function. This permits the designer to optimize gas turbines 
either for weight or for losses depending on the specific 
purpose of the machine. By giving different weights to the 
individual objective functions, a parametric study can be 
conducted to find the relative effects of the two objectives on 
the optimum design. 

The analysis used in this work accounts for the 
Reynolds n-umber effects and makes use of free vortex theory to 
include three dimensional effects in turbines. The various 
behaviour constraints which arise due to fluid flow for the 
efficient operation of the stage are included. Rrom mechanical 
design point of view, constraints are placed on the deflection, 
stress and fundamental natirral frequency of vibration of the 
blade . 

A finite element procedure is developed for finding 
the stresses, deflection and natiiral frequencies of vibration 
of turbine blades. This procedure takes into accoimt the 
considerations like taper, twist, asymmetry, rotation, rotary 
inertia and shear deformation of the blade along with the gas 
loading. Although the idealization of gas turbine blades can 
also be made by using more complex elements like flat or 
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curved triangular plate elements, the beam element has been 
used in the present work as it gives reasonably good results 
with lesser number of degrees of Ireedom and requires lesser 
computer time. 

In the optimization of complex problems using finite 

element method, a designer is generally confronted with two 

problems, namely, the computer storage and the computer time. 

To economize computer time, the eigen value problem has been 

solved by using one of the most efficient solution techniques 

8 1 

developed by Bathe and Y/ilson for large structural systems. 

In Bathe and Wilson technique, the Rayleigh-Ritz sub-space 
iteration algorithm, which solves the eigen value problem 
directly without transformation to the standard form, is used. 
Before solving the eigen value problem, the deflection and 
stresses due to gas forces have been calculated to make use of 
the available stiffness matrix. In this work, the Cholesky 
decomposition of symmetric matrices, storing only the upper 
triangular matrix, is used for solving the equilibrium equations. 

In order to avoid storage of the necessary air 
properties in tabular form, the air properties have been 
computerised with polynomial approximations using least squares 
technique. 

The constrained optimum design problem has been cast 
as a nonlinear mathematical programming problem. The interior 
penalty function method, with a variable metric unconstrained 
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minimization technique, is used to solve the optimum design 
prohlem. The minimizing step lengths in the unconstrained 
minimization are determined by cubic interpolation method. 

The computation of the gradient and the slope of the functions 
has been carried out by using backward difference method. 

1 .3 Organization of the Thesis 

The thesis is organized into eight chapters and five 
appendices. After the introductory chapter, the design 
philosophy, objective function, design constraints and design 
variables of the problem are stated in Chapter 2. The problem 
of optimum design of axial flow gas turbine stage is also 
formulated in this chapter. 

The evaluation of the objective function, which 
includes derivation of expressions for the weight of stage and 
the losses in stage is considered in Chapter 3 . The iterative 
procedure involved in finding the efficiency (one minus losses) 
from the known data and the computation of the geometrical 
properties of airfoil section are also discussed in this 
chapter. 

The development of finite element procedure for 
analysing turbine blades with consideration of taper, twist, 
rotation, rotary inertia and shear deformation is given in 
Chapter 4. The method of evaluating the constraints of the 
optimization problem is discussed in Chapter 5. A procedure for 
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converting airfoil cross section into an equivalent rectangular 
section is also presented in this chapter. 

Chapter 6 deals with the optimization algorithm used 
in the present work. The reasons for using the penalty 
function formulation, with a variahle metric unconstrained 
minimization method, along with the computational details are 
given in this chapter. 

Several illustrative examples are presented in 
Chapter 7 to show the effectiveness of the present approach 
for the design of axial flow turbines. The results of sensi- 
tivity analysis are also presented in this chapter. The 
conclusions drawn from the present work and the scope for 
further extension of the work are stated in Chapter 8. 

The polynomial equations used in the calculation of 
profile losses are given in Appendix A. The polynomial 
equations for evaluating the various air properties ^e stated 
in Appendix B. Appendix C gives the formulation of various 
matrices required in evaluating element stiffness and mass 
matrices of Chapter 4. In Appendix I), the elements of stiffness 
and mass matrices of a tapered beam element are presented. 
Appendix E contains a description and other details of the 
computer programme developed for the optimum design of gas 


turbine rotor stage. 
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CHAPTER 2 

POmiULATIOlf OE THE OPTIMIZATION PROBLEM 

Vfhen a means for predicting the behaTior of any design 
within a particular design concept is available, limitations on 
the performance and other external constraints on the design are 
stated and an acceptance criterion is established, it is possible 
to cast the design mod'’’f ication problem in the form of a mathe- 
matical programming problem. 

2.1 Design Philosophy 

Any general design problem can be formulated and solved 
accoiding to either deterministic or probabilistic design philo- 
sophy. If all the quantities affecting the design problem are 
deterministic, the design problem can be solved according to the 
deterministic design philosophy. On the other hand, if some of 
the design parameters are random in nature , the design problem 
has to be solved according to the probabilistic design philosophy. 

In the deterministic design philosophy, a general 
mathematical programming problem can be stated as follows; 

Minimize a multivariable function f(X) where X is a 
n-dimensional vector consisting of X., j = 1, 2, ..., n, subjected 
to the given constraints g^(X){£ ,=,i^}0, i=:1,2, ...,m. 

The function f (X) is called the objective function and the 
vector X is termed as design vector. 
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In the probabilistic design philosophy, a general 
mathematical programming , problem can be stated as follows; 

Minimize the multivariable function f(X, Y) subjected 
to the probabilistic constraints P[ g. (X, Y){^, } 0] ^ s 

j = 1, 2, m, where X is the vector of design variables, and 

Y is the vector of other parameters affecting the design problem- 

-y -y 

Here the components of X and Y are assumed to be random variables, 

f represents the mean value of the objective function and 

P[ g. (X, Y){<^,=,^}0] >p. denotes that the j constraint 

has to be satisfied with a probability of greater than or equal 

to some specified quantity, p., j = 1, 2, ..., m, where 

J 

® 1. Pi i * 

In the present work, only the deterministic design 
philosophy is used for the optimization of axial flow gas 
turbine stage . 

2.2 Objective Punction 

A design problem usually has several solutions which 
may satisfy the specified functional requirements adequately. 

The objective function in a general optimization problem 
represents a basis for choice between alternate acceptable 
designs. In most of the practical design problems the mini- 
mization of weight, cost, volume or losses, or maximization of 
profit, rigidity or efficiency is taken as the objective. In 
the case of gas turbines used in aerospace applications the 
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minimization of weight is one of the most important criteria 
while in. 1 iie case of gas turbines 'used in stationary power 
plants, the maximization of efficiency represents a more useful 
criterion. In some cases a mixed objective function represen- 
ting a linear combination of weight and efficiency will be a 
more useful objective. In this work the mixed objective 
function is used so that the optimum design of aerospace turbines 
or industrial turbines can be obtained from the same computer 
program going ‘ proper weightages to weight and efficiency in 
the objective function. 

2.3 Design Constraints 

In a realistic gas turbine design generally the 
following requirements are to be met from aerodynamic, vibra- 
tional and strength considerations: 

The rotational velocity of the rotor should be within 
some upper and lower bounds . 

The aspect ratio (height/chord) of rotor and nozzle 
blades should be within some specified upper and lower 
bounds . 

The pitch-chord ratio of the rotor and nozzle blades 
should be within the specified upper and lower bounds. 

The relative gas velocity angles at inlet and outlet 
of rotor blades ^2 ^3 should be within some 

specified upper and lower bounds . 


(i) 

(ii) 

(iii) 

(iv) 
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(v) The axial velocity of flow should be within some upper 
and lower bounds, 

(vi) The actual pressure ratio across the nozzle blades 
should be below the critical pressure ratio. 

(vii) The Mach number at the exit from the stage should be 
less than a specified value. 

(viii) The included angle of divergence of the turbine annulus 
walls should not exceed the specific upper limit. 

(ix) The flow coefficient { 0 ) and the stage temperature drop 
coefficient ( 4 > ) should lie within some upper and lower 
bounds . 

(x) The degree of reaction at mean radius should be within 

some upper and lower bounds while the degree of reaction 
at the root of rotor blades should be a non-negative 
quantity . 

(xi) The fundamental natural frequency of the blades should 
be away from the forcing frequency of the blade in 
order to avoid resonance. 

(xii) The stresses developed at the root of the rotor blade 
should be less than the pemissible value. 

(xiii) The tip deflection of the rotor blade should be less 
than some specific value. 
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2.4 Design Variables 


Once the objective function and the design constraints 
are specified the mathematical programming problem can be 
stated as soon as the design variables are identified. Dor the 
design of a gas turbine stage, the following parameters are 
taken as design variables: 


X 


1 


X, 


X. 


X. 


X 


6 


X 


7 


"8 


= Mean diameter of the rotor, d 

= Ratio of the chord of rotor blade to mean diameter 


'R 


’ d 


= Ratio of the chord of nozzle blade to mean diameter. 


'H 


X. = 


Ratio of spacing to diameter at mean radius of the 

s. 


nozzle blades, 




= Ratio of spacing to diameter at mean radius of the 

®R 

rotor blades, ^ 

= Relative angle of the velocity triangle at the inlet of 
the rotor blade at mean radius, $2 
= Relative angle of the velocity triangle at the exit of 
the rotor blade at mean radius, 

= Axial velocity of flow across the stage, C . 

Thus the vector of design variables X becomes 


X 




X 

X, 

£ 

x' 

X, 

X, 


1 


4 ^ 


7 

^8 


V. 


Oj,/a 

Sjj/d 

Sj/d 


(2.1) 
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It can be noted that the nozzle angle can be taken 
as a design -variable in place of the axial velocity as they 
are related by the equation 

C = V ^ y- (2.2) 

a tan ctg “ ^2 

where U is the tangential velocity of the rotor at mean radius. 


2.5 Optimization Problem 

The optimization problem can now be stated as follows; 

Minimize 


f(i) = K^d - 

+ ^2 "m^ ■ 

2(1 

(a-hjj)2 

+ B ) 






( 2 . 3 ) 

subject to 





Si = 

u(i) 

U - 1 .0 5 

0 


( 2 . 4 ) 

sp = 

- 1 .0 < 

0 


( 2 . 5 ) 

g3 = 

S (1) °TJ 

(_S.) _ (_R) 

< 

0 

( 2 . 6 ) 

§4 = 

% % 

(_R) _ (_I) 

< 

0 

( 2 . 7 ) 

S 5 = 

St 

(_1) _ (J) 

< 

0 

( 2 . 8 ) 

S6 = 

c (u) 

(^) _ 

S ^ ^d ^ 

< 

0 

( 2 . 9 ) 
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S? = 

S-^T (l) S-nJ 

(^) - (5^) < 

0 

(2.10) 

% = 

^TCT ^TT 

{^) - (#) i 

0 

(2.11 ) 

gg = 

S„ (1) Sp 

(3^) - (-6-) i 

°E R 

0 

(2.12) 

§10 

S-p Sp (u.) 

(JS) - (^) < 

R R 

0 

(2.13) 

§11 

4^^ - ^2 - ° 


(2.14) 

§12 

&2 - ^2"^^ 1 0 


(2.15) 

§13 

- K 10 

3 3 


(2.16) 

§14 = 

- 6(«) i 0 


(2.17) 

§15 "" 

0(1) 

^ - 1 .0 <0 


(2.18) 

§16 

C 

■ — ^ _ 1 .0 < 0 
q(u) 


(2.19) 





§17 ^ 

Po1 . £01 ,0 

P2 Pc ”■ 


(2.20) 

§18 = 

0 

vl 

;:$ K\ 

1 

K> 


(2.21 ) 

§19 = 

a - < 0 

C C 


(2.22) 

§20 

0 ^^^ - 0 lO 


(2.23) 
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§21 

§22 

§23 

§24 

§25 

§26 

§27 

§28 

§29 

§50 


0 _ < 0 


- 4i < 0 


4 . - < 0 

- E < 0 


R - R 




a _ ot < 


2 - 


R^^^, - R . < 0 

root root — 


- 1 .0 <0 


- < 


- 1.0 < 0 


(2.24) 

(2.25) 

( 2 . 26 ) 

(2.27) 

(2.28) 

(2.29) 

( 2 . 30 ) 

(2.31) 

(2.32) 

(2.33) 


where 

f(X) = mixed objective function 

= weightage given to the minimization of losses in 
the objective function 
Hg = stage efficiency 

= weightage given to the minimization of weight of 
stage in the objective function 
p = mass density of disc and blade material 

3 = chord taper ratio 
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= chord of rotor blade at mean radius 
h„ = mean height of rotor blades 

It 

ru = number of rotor blades 

K 

= cross-sectional area of rotor blade at the mean 
radius 

Ujj = number of nozzle blades 
h|jj- = mean height of nozzle blades 

Ajj = cross-sectional area of nozzle blade at the mean 


®R 

Pol 

P2 

Pc 

]VU 

^5 


a 


c 


R 

2 

R 4. 

root 


radius 

chord of nozzle blades at mean radius 
spacing of nozzle blades at mean radius 
spacing of rotor blades at mean radius 
stagnation pressure at the inlet of the stage 
static pressure at the inlet of the rotor blades 
critical pressure 

Mach number at the exit of the rotor blade 
calculated from the absolute velocity 
included angle of divergence of flared turbine 
annulus 

degree of reaction at mean radius 
nozzle blade angle at the outlet 

degree of reaction at the root of the rotor blade 
rotational speed of the rotor (revolutions per 
second) 
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£i5(l) = first natural frequency of vibration of rotor 

blade (cycles per second) 

0 ^ ^ = maximum stress at the root 

= tip deflection of the rotor blade 
1 = superscript denoting the lower bound 

u = superscript denoting the upper bound 

The quantities Pg , , 0, » R, ^2 ’ \oot’ ^t 

and stresses are dependent on the design variables and hence 
the constraint . equations (2.20) to (2.33) ame called the behavior 
constraints. The relation between the behavior quantities and 
the design variables cannot be expressed directly in closed 
form. However, for any given design vector X, the behavior 
quantities can be evaluated numerically upto any desired 
accuracy. The analysis procedure to be adopted for determining 
the behavior quantities is considered in Chapters 35 4 and 5. 

Equations (2.4) to (2.19) represent the geometrical 
or side constraints, which impose limits on the size of the 
design variables. 

It can be seen that the objective function of equation 
(2.3) is a non-linear function of the design variables, and the 
side constraints of equations (2.4) to (2.19) are linear while 
the behavior constraints are non-linear. Hence the mathematical 
programming problem formulated above is a non-linear program- 
ming problem. 
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CHAPTER 3 


EVALUATION OE OBJECTIVE FUNCTION 
OF AXIAL FLOW GAS TURBINE STAGE 


The essential purpose of a turbine is to obtain work 
output from a gas which drops in pressure and hence in tempera- 
ture while passing through the machine . There are two basic 
types of gas turbines, namely, radial flow and axial flow 
turbines. The radial Turbine is similar in appearance to a 
centrifugal compressor, but with a flow in the inward direction 
and nozzle vanes replacing the diffuser vanes. The radial 
turbine has been used in some small gas turbines where compact- 
ness is more important than performance. In an axial flow 
turbine, the gas flow has components of velocity parallel to the 
axial and tangential directions of the machine but has little 
radial velocity. The range of applications of axial flow 
turbines is very wide. There is no alternative to the axial 
flow turbine for the designer who seeks to have high flow per 
unit frontal area and high efficiency coupled with a reasonably 
high enthalpy drop per stage. In axial flow turbine designs, we 
want the machine to have a higher efficiency with a lesser 
weight. That is why the losses of the turbine stage and its 
weight have been taken as objectives in the present woik. The 
procedure for evaluating these objective functions is discussed 
in some detail in this chapter. 
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3.1 Two Dimensional Fluid Flow Idealization in a Gas Turbine 
Stage 

Figure 3.1(a) shows a gas turbine stage consisting of 
two rows of blades, one stationary and the other rotating. In 
the 'stator' or 'nozzle' row, the tangential -velocity of the 
fluid is increased in the direction of rotation and the fluid is 
accelerated at the expense of a pressure drop. The tangential 
velocity in the direction of rotation is reduced across the 
'rotor' row and tangential forces are exerted by the fluid on 
the rotor blades which produce a torque on the output shaft. 

Thus the absolute kinetic energy of the fluid is reduced across 
the rotor. There may be some acceleration of the fluid relative 
to the moving blades together with some drop in pressure and 
temperature. A multistage turbine is made up of several stages, 
each stage consisting of first a nozzle row and then a rotor 
row . 

For machines having high hub to tip ratio, there cannot 
be large radial components of velocity between the annular walls 
and the flow conditions are not very much different at different 
radii. For such machines the two dimensional analysis is 
sufficiently accurate. In two dimensional analysis, the flow 
velocity triangles are drawn for the conditions of mean radius, 
but these triangles are assumed to be valid for other radial 
sections as well. 
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3.2 Velocity Triangles 

Figure 3 . 1 (l) shows the Telocity triangle for an axial 
flow turbine stage alongwith the nomenclature employed. The gas 
enters the row of nozzle blades with a static pressure p^ » 
temperature T^ and velocity , expands to P2 ? ’^2 leaves with 
an increased velocity Cg at an angle 02* The rotor blade inlet 
angle will be chosen to suit the direction B2 velocity 

V2 relative to the blade at inlet. Bg and are found by 
vectorial subtraction of the blade speed U from the absolute 
velocity C2. After being deflected, and usually further expanded, 
in the rotor blade passages, the gas leaves at p^, T^ with 
relative velocity V^ at an angle 3 ^. The vectorial addition of 
U and V^ yields the magnitude and direction of the gas velocity 
at exit from the stage, and o^. The angle is known as 
the swirl angle . 

In a single stage turbine, will be axial, i.e. 
a. = 0 and C. = C where C is the axial velocity of flow. If 
on the other hand, the stage is typical of many in a multi stage 
turbine, and will prob^ly be equal to and so that 
the same blade shapes can be used in successive stages. The 
quantity (C + C ) represents the change in whirl (or tangen- 

W2 w^ 

tial) component of momentum per unit mass flow which produces 
the useful torque. 
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3 ..3 Temperature -Entropy (T - S) Diagram for a Stage 

Figure 3*2 represents the temperature-entropy diagram 
for a stage. The full and chain dotted lines connect the 
stagnation and static states respectively. The stagnation temper- 
ature T ^2 equal to T^^ because no work is done in the nozzles; 
and the short horizontal portion of the full line represents the 
stagnation pressure drop, (p^^ - friction in nozzles. 

The losses are of course exaggerated in this figure. When obta- 
ining the temperature equivalent of the velocity of the gas 
leaving the nozzle row, we may say that ideally the gas would be 
expanded from T^^ to T^ but due to friction the temperature at 
the nozzle exit will be T^ which is somewhat higher than T^. 

The expansion of the gas in the moving blade passage reduces its 
pressure to p^ . Isentropic expansion in the whole stage would 
result in a final temperature T^ , and in the rotor blade passages 
alone T^ . The expansion with friction leads to a final temper- 
ature T^. As no work is done by the gas relative to the blades, 
the stagnation temperature relative to the blade at point 3 , 

Tq^ rel’ equal to T^^ stagnation pressures 

relative to the blade at inlet and outlet are represented by 
Po2 rel Po3 rel '^espeotlvely. 



32 


3.4 Losses and Efficiency 

Horlock®^ has discussed about losses and efficie’^'^y of 
full admission turbines and Yahya ^ has discussed about losses 
in axial flow turbines with partial admission. Here the infor- 
mation regarding losses and efficiency of axial flow turbines 
with full admission is summarized. 

3.4.1 Blade losses due to Blow of Fluid in an ixial Flow 
Turbine Stage 

An overall blade loss coefficient Y (or ) must account 
for the following sources of friction loss: 

(a) Profile loss - associated with the boundary layer growth 
over the blade profile (including separation loss under 
adverse conditions of extreme angles of incidence or high 
inlet Mach number) . 

(b) Annulus loss - associated with the boundary layer growth on 
the inner and outer walls of the annulus shown in Figure 
3.3(a). 

(c) Secondary flow loss - arising from the secondary flows which 
are always present when a wall boundary layer is turned 
through an angle by an adjacent curved surface. This loss 
is shown in Figure 3.3(b). 

(d ) Tip clearance loss - near the rotor blade tip the gas does 
not follow the intended path, fails to contribute its quota 
of work output and interacts with the outer wall boundary 
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layer. The radial tip clearance and side clearance for 
shiouded hlades are shown in figures 3.3(c) and 3.3(d) 
respectively. 

The profile loss coefficient is measured directly in 

the cascade tests. The annulus and secondary flow losses cannot 

be easily separated , and they are accounted for by a secondary 

loss coefficient Y . The tip clearance loss coefficient, which 

s 

normally arises only for rotor blades, will be denoted by Y^^. 
Thus the total loss coefficient Y comprises of the accurately 
measured two dimensional loss Y^, plus the three dimensional 
loss (Y_ + Y, ) which must be deduced from the turbine stage test 

S 

results. 


3.4.2 Definition of the Loss Coefficients 

O / 

Benson ^ reviewed the methods for assessing the loss 

coefficients in radial gas turbines while the various loss 

coefficients used in axial flow compressors and gas turbine 

designs have been defined, compared and discussed in detail by 
83 

Brown . Referring to Figure 3.2, the most common loss coeffi- 
cient for the nozzle blades msy be defined either by 


(o|/20p) 


or 




PqI - Po2 
Pol - P2 


(3.1) 


Both X and Y express the proportion of the leaving energy which 
is degraded by friction. Although Yjj can be measured relatively 
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easily in cascade tests, can be used more easily in designs. 

It can b3 shown that and Yjjj. are not very different numerically. 
Similarly the loss coefficient for the rotor blades is defined as 


I- 


R 


T" 
3 


(V^/2Cp) 


or 


= 


^o2 rel ~ ^o3 rel 
^o3 rel ~ ^3 


(3.2) 


3.4.3 Definition of Efficiency 

The isentropjc stage efficiency or the total -to-total 

stage efficiency ri_ is defined as 

s 


T - T 


(3.3) 


From Figure 3.1(c) JLt can be seen that 


^o3 " ^<^3 " ^^3 ■ ^3^ = ^®3 ■ ^3^ + ^^3 - ^3^ 


(3.4) 


But (T^/T^) = (Tg/T^) because both of them are equal to 
( Y -1 ) 

(~) . Rearranging and subtracting one from both sides of 

P2 

equation (3.4), one gets 


T” — T ^ 

3 3 

T’ 

3 


T - T ' 

12 

T ' 

2 


or 


(T^ - T^) 


T 

(Tq - ^ (3.5) 


Hence 
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1 


1 + [{(T^ - ^3) + (^2 




(3 .6 ) 


1 




(3.7) 


Alternatively, by substituting V_ = C sec 3 C» = C sec « , 

^ B. ^ £L Ei id 

and 


- '^03^ 


U C (tan 3, + tan B^) 

3. j d. 


U C , tan 8_ + tan a„ 
aj 3 2 


(SL) 

‘ f 
a J 


(3.8) 


equation (3.7) can be written in the form 


"s = 


1 


T 

^3 (j^) 3jj sec^ 

tan + tan Og - (^) 

8. 


1.1^ 

‘ ^ 2 u 


(3.9) 


Because Y - X , the loss coefficients Yj^ and may replace X 
and X^ in equations (3.7) and (3.9) if desired. 

A second definition, the total -to-static efficiency, 
n^g, for a stage is given by 


R 


TS 


5^01 - 


( 3 . 10 ) 
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The isentropic stage effic lency should generally be 
used in multi stage turbines where the exhaust Telocity from a 
stage is not lost, 

3.5 Evaluation of the Efficiency of a Gas Turbine Stage 

In Chapter 2 we identified eight design variables at 
mean radius as; diameter of the rotor, d; chord to diameter 

^T) 

ratio of the rotor blades, * 3 — I chord to diameter ratio of the 
nozzle blades, j--, spacing to diameter ratio of nozzle blades, 

. . . Sp 

^ ; spacing to diameter ratio of rotor blades, — ; gas angle for 
rotor blade at inlet, angle for rotor blade at outlet, 

and the axial velocity of flow across the stage, C^. In 
optimization procedure, the efficiency of the turbine stage is 
to be evaluated at various trial values of these design variables. 

In the design of axial flow gas turbine stage it is 
assumed that the stagnation pressure at inlet to stage ) » 
stagnation temperature at inlet to stage mass flow rate 

across the stage (m) and the speed of the rotor (N in revolutions 
per second) are the preassigned parameters. 

Further we assume that the properties of air like 
specific heat at constat pressure, c , gas constant, E, visco- 

Jr 

sity of gas, p, and ratio of specific heat y are known or can 
be calculated from suitable formulas. 
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In the following suhsectio'ns the equations necessary 
for evaluating the isentropic stage efficiency or total-to- 
total stage efficiency n_ are developed sequentially. It is 

o 

essential to assume a trial value in the beginning to calculate 
a more accurate value of the efficiency ri^. Thus the process 

is iterative and has to be continued until the values of n in 

' s 

two consecutive iterations are sufficiently close to each other. 

3 .5.1 Calculation of the Parameters of Telocity Triangle and 
Blade Heights 

The tangential velocity of the rotor, U, can be 
calculated as 

U = TTdN (3.11) 

With the known values of U, G , 3^ and 6 , the velocity tri- 

angles (Pigure 3.1(b)) at inlet and outlet of the rotor can 
be drawn and their parameters can be calculated using the 
following equations once the value of C is assumed constant 
across the stage. 

C 

0 ^ ( 3 . 12 ) 


R - ^ (tan 3^ - 

tan $ 2 ^ 

(3.13) 

tan = tan ^2 


(3.14) 

tan = tan 3, 

3 3 

- F 

(3.15) 
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C 

a 

cos <X2 
C 

a 

cos a- 
3 

C 

a 

cos ^2 
C 

a 

cos 

3 

°2 “2 

sin a_ 

J D 


(3.16) 

(3.17) 

(3.18) 
(3 .19) 
( 3 . 20 ) 
( 3 . 21 ) 


where 0 ani R represent the flow coefficient and the degree of 
reaction respectively. 

The stagnation temperature drop in the stage 

and the blade loading coefficient or the temperature drop 
coefficient can be calculated as 


A T 


TJ G (tan 3 ^ + tan 0 „) 
a d j 


os 


( 3 . 22 ) 


and 


Tf. = 


2c AT 
p os 


(3.23) 


The temperature equivalent of the outlet velocity is -5— and 

^ P 

hence 


m m _ 

^2 o2 2c, 


= T 


o1 2c 


P 


(3.24) 
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Assuming a guess value for the loss coefficient of the nozzle 
■blades, equation (3.1) gives 


c; 


T ’ 

2 


T - 
2 


N 2c. 


(3.25) 


The pressure P 2 can "be found from the isentropic relation as 

y 




( 3 . 26 ) 


By ignoring the effect of friction, the critical pressure ratio 
can be calculated as 



Y + 1 


Y 

(y-i"J 

) 


(3.27) 


where p is the critical pressure. The actual pressure ratio 
Pol 

— ^ is kept well below the critical pressure ratio for avoiding 
P2 

choking of the nozzle. The density Pg? annulus area and 
the throat area of the nozzle A 2 JJ are calculated from the 
relations: 


and 


2 “ — ’ 

R T 2 


m 


"a "a ’ 


m 


"2R 


2 ^2 


Co “ ^ 


COS Ot, 


(3.28) 

(3.29) 


(3.30) 
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Y/e now assume that the C. is axial i.e. C. = C and 

1 1 a 

hence Gc- = 0. The temperature equi'valent of the inlet kinetic 
' ? 

energy is and hence the annulus area required in plane 1 

P 

can be calculated as follows: 


and 


T. 


Pi 

Pi 



'o1 



Y 


rn C y y 

Ol 


p- 


R T. 




m 



(3.31 ) 


(3.32) 

(3.33) 


(3.34) 


Similarly at the outlet of the stage, we have 


T 


o3 


= T . - i T , 
ol os ’ 


T. 


T _ ^ 

» 


at 


Y 

Ty=TT 


'o3 




T TyITT 
Po3^T ^ 


o3 


(3.35) 

(3.36) 

(3.37) 

(3.38) 
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^3 



and 


(3.39) 


A. 


m 


a 


(3.40) 


The rotor blade height h^^ and the nozzle blade height 
hjj can be calculated from the blade heights h^ , h^ and h^ at 
planes 1 , 2 and 3 using the relations 


^1 

A^N 

U 

(3.41 ) 

^2 

1 

1 

II 

(3.42) 


II 

(3.43) 


^2 ■*' ^3 

2 

(3.44) 


h^ + 

2 

(3.45) 


The turbine annulus is flared as shown in Figure 3.4. 

It is desired that the included angle of diTergence of the wall 

a should not exceed a certain limiting value. If the blade 
c 

breadth is approximately taken as the chord of the blade , 
is related to the other parameters as 



(3.46) 
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It is necessaiy to check the Mach number at the exit 

of the stage, Mp , because if this is too high, the friction 

3 

losses may become unduly large . The Mach number at the exit 
is computed as 


C, 


M. 


( Y RT^)^ 


The loss coefficient can be calculated as 




T_ - T" 
_2 2 

V^/2c 

3^ P 


v/here 


T^ 


Po 


( 1 - 1 ) 

Y 


(3.47) 


(3.48) 


(3.49) 


3.5.2 Calculation of Gas Angles Along the Radius by Considering 
Three-Dimensional Effects 

In actual practice, the shape of the velocity triangle 
varies from root to tip of the blade since the blade speed U 
increases along with the radius.. Another reason is that the 
whirl component in the flow at outlet from the nozzles causes 
the static pressure and tmperature to vary across the annulus. 
With a uniform pressure at inlet or with a small variation in 
pressure (due to smaller whirl component), it is clear that 
the pressure drop across the nozzle will vary giving rise to 
a corresponding variation in efflux velocity C2. Twisted 
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blading designed to take account of the changing gas angles is 
called vortex blading. 

It can be shown that if xhe elements of the fluid are 

to be in radial equilibrium (— ^ = — ) , an increase in static 
pressure from root to tip is necessary whenever there is a 
whirl component of velocity. The gas turbine designer cannot 
talk of impulse or 50 percent reaction stages as the rotor 
stage pressure or temperature drop increases from root to tip 
and therefore the reaction will also increase from the root to 
the tip. 


The design method used in the present analysis, which 
takes into account these three dimensional effects, is called 
the radial equilibrium free vortex design. The assumptions of 
this analysis are; 

(i) The stagnation enthalpy h is constant over the annulus 
dh ° 

“ = 0 ). 


(i.e . 


dr 


(ii) The axial velocity is constant over the annulus. 

(iii) The whirl velocity is inversely proportional to the radius. 

By using the subscript r to represent a quantity at 
any general radius r, the following relations can be derived 

for calculating the angles “2r’ °'3r’ ^2r Sr 
general radius ; 


tan a 


2r 


(^) tan « 
^r 2 


( 3 . 50 ) 



44 


oail a. 


tan 3 


(f-) ■tan a 
°-r 3 ^ 




2r 


tan B 


(|-) tan + (^) ^ 


rN U 


3r 


(3.51 ) 

(3.52) 

(3..53) 


‘r 3 


3 a 


These general relations can be specialized to obtain 
the blade angles at the root and at the tip of the stage as; 


d 


i^aii <^2 root ^d - h^ 


■) tan a. 


(3.54) 


xan a 


3 root 


tan a 


2 root 


tan a 


3 root 


“3 


/ a ^ ^ 4 ” ^ 2 ^ ti 

“2 - <-X — ) — 

Si 

3 a 


(3.55) 

(3.56) 

(3.57) 


tan a, 


tan a 


tan a 


2 tip 


3 tip 


2 tip 


“ 2 


“3 


, d N , ^ ^2 3 IL 

^d + h^^ “2 " ^ d ^ C 

2 a 


tan B 


3 tip M + h. 


u ) tan + (“ 


d + h. 




(3.58) 

(3.59) 

(3.60) 

(3.61) 


a 
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It is desirable always to have a positive value of the 
degree of reaction at the root of the blade, ^poot ’ 
be calculated from the following relation: 


root 


{ tan 


3 root 


- tan 0 


2root 


(3.62) 


3.5.3 Blade Angles and Blade Profile 

It is important to note that the velocity triangles 
yield the gas angles but not the blade angles. We have shown 
how to establish the gas angles at all radii and the blade 
height. We have also taken some trial value of spacing and 
chord of the blade. The next step is to choose the stator and 
the rotor blade shapes which will accept the gas incidence upon 
the leading edge and deflect it through the required angle with 
minimum loss. For calculating blade angles we should know the 
value of deviation. 

86 

Carter and Hughes have calculated the deviation (5) 
in potential flow through a variety of compressor and turbine 
cascades. The cascade geometry is shown in Figure 3.5. It 
has been suggested that the un-stalled low speed deviation for 
the turbine cascades may be written as 


5 = m 0' v'f (3.63) 

w 

where m is a function of the stagger a as indicated in 

0 

Figure 3.6(a) and 6' is the blade camber angle. The angle 
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9-’ is related to the blade aiigles (a^ and a^) of the cascade 
and the gas angles of the cascade ( and at inlet and 
outlet as 


9' = + i) + ( + 6 ) = ^ + ± + S (3.64 ) 

where i is the angle of incidence and e is the gas deflection 
which is equal to (a^ + a^). At design conditions, the angle of 
incidence is assumed to be zero and hence the equation of 
deviation becomes 


5 = + ct^ + 6 ) /| (3.65) 

The solution of equation (3.65) gives the deviation as 


m( + a^) 
(1 - m /|) 


(3.66) 


In the present analysis, a constant value of m = 0.1 
has been assumed which is reasonably good for the blade having 
parabolic arc camber line. Thus we can calculate the blade 
angles at the outlet of the nozzle stage and the outlet of 

the rotor blade from the relations given below. The inlet 
stator and rotor blade angles will then correspond to and 
respectively as incidence is asstimed to be zero. 


06 2 


s = 


Clr 


+ 


m( 



(1 


- /®isr 

m V 


) 


«2 


(3.67) 
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m( Cg + 

0 ^ = + ( 3 . 68 ) 

0 -i/^) 

The blade angles at the root ani the tip for free 
vortex blading can be calculated by replacing and ct^ of 
equations (3.54) to (3.6l) by a ^ and of equations (3.67) and 
( 3 . 68 ). 

The next step is to fix the blade profiles. Two 

series of profiles have been specially designed for turbines 

8T 

by Dunavant and Erwin . The details, of the two series are 

given in Table 3.1 and in Eigure 3.7. The primary series A^K,^ 

(Figure 3.7(a)) is for reaction blades in which there is 

acceleration through the cascades. The shape of the camber 

line gives a rapid turning in the f oiward part of the blade , 

where the Mach numbers are low. The secondary series 

(Eigure 3.7(b)) is for conditions approaching zero reaction. 

The loading is distributed more evenly along the chord, since 

there is less change in the average Mach niimber through the 

channel. The camber line coordinates and the blade thickness 

y, have been given as a percentage of the chord in Table 3.1 
^ dy 

for which = 0.8574, leading edge radius = 4.407^ of 

iy ^^ 0 ( 0 . 5 ) 

Chord; = -0.1602, trailing edge radius = 1 , 000 % of 

^^c(95) 

chord for primary blade and for secondary blade 
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dye 

dx (0.5} 


0.5657, 


_ 

dx^(95) " 


- 0 . 2017 , 


leading edge radius = 3.300% of chord j and 
trailing edge radius = 1.000% of chord. 


The abo-ve -mentioned profiles and are 

suitable for a particular set of inlet and outlet angles. If 
the gas angles are changed, we will have to select another set 
of profiles whose mean line coordinate should match with the 
new gas angles . In this way we will have to store a huge amount 
of data in an automated design procedure which may not be 
practically feasible . This problem can be solved if we find 
the mean line of the profile by satisfying certain desired 
conditions and then superimpose the thickness distribution of 
either the primary or the secondary profiles from Table 3.1 
depending on the requirements. 

In Figure 3.8, the line 1-2 represents the chord of 

the blade and the curve 1-m-2 the camber where the maximum 

camber with respect to the chord occurs at the point m. The 

angle between the chord line (x-axis) and the axial direction 

(X-axis) is given by the stagger angle a . We assume the loca- 

s 

tion of the point m depending on the type of blading (reaction 
or impulse) and hence the coordinates of m (namely, Xj^^ and y^) 
are known. It is required that the gas should enter at point 1 
(when y = y^ ) at an angle a' and it should leave at point 2 
(when y = y^) at an angle og* Hence the following conditions 
are to be satisfied by the camber line of Figure 3.8; 
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At 

X 

= : 

y = y-i 

and 

il 

a 

1 


At 

X = X2 

= : 

y = y2 

and 

^ _ 
dx 

"2 

(3.69) 

At 

m 

0 

y = ^m 

and 

ll 

0 



Thus there are six conditions to be met by the camber 
line and hence a polynomial of degree five can be fitted 
exactly. Once the equation of the mean line is obtained, we 
can superimpose the thickness distribution from Table 3.1 to 
obtain the required profile of the blade . The upper and lower 
profile curves of the blade section have been interpolated as 
fourth order polynomials from 23 data points using least 
squares method of curve fitting. These boundary profiles have 
been expressed as y^^^(x) and y^^\x) respectively in Fig. 3.8. 

3 . 5.4 Estimation of Design Point Performance of the Stage 

Once the major parameters of the stage are' decided, 

the design point performance of the stage can be calculated. 

88 

The method used is due to Ainby aM Mathieson which estimates 
the performance on flow conditions at the mean diameter of the 
annulus. Reference 88 describes how to calculate the perfor- 
mance of a turbine over a range of operating conditions, but 
we shall be concerned here only to find efficiency at the 
design point. A start is made using the two correlations for 
profile loss coefficient Yp obtained from cascade data (shown 
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in Figures 3.6(b) and 3.6(c)). These refer to nozzle-type 
blades ( 02 = impulse type blades ( 02 “ ^3^ conven- 

tional profiles having thickness to chord ratio (— ) of 0.2 and 

t ^ 

a trailing edge thickness to pitch ratio, (") ? of 0.02. The 

o 

rotar blade notation is used in Figure 3.6 and the flow relative 
to any blade row is considered. When the nozzle row is being 
considered, becomes and becomes The values of 

in Figure 3.6 refer to the blades operating at zero incidence, 
i.e. when the gas inlet angle 02 also the blade inlet angle. 
The incidence is assumed to be zero and the thickness to chord 
ratio of the rotor and nozzle blades is also assumed as 0.2. 

All the curves of Figures 3.6(b) and 3.6(c) have been comput- 
erised in the form of polynomials. The polynomial equations 
of the various curves are given in Appendix A. The method of 
estimation of the design point performance can be stated by 
the following steps. 

Step (i): Estimate the profile losses of rotor and nozzle 

blades as follows; 


P R 


02.2 


^p(0o= 0) " 


■p( 62 = 0)^ 


tTj/cTj ^2^ ^3 
. R^ E x 

0»2 ^ 


(3.70) 
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,( a^ = 0)}3 ^ 0.2 


%/°U- 2 


) 


(5.71) 


s 


Step (ii) ; The secondary and tip clearance loss data for Y 
and have been correlated using the concepts of lift and drag 
coefficients. By defining the mean angles a ard 3^ as 


tan a 


m 


(tan ® 2 ~ ^ 

2 


( 3 . 72 ) 


and 


tan 0 


(tan 3^ - tan 3^) 


m 


( 3 . 73 ) 


the lift coefficients for nozzle and rotor are calculated from 
the relations 


®-RT 

(Gt ) = 2(— )(tan a + tan a ) cos 

^ B Cu 1 2 

and 


a 


m 


(3.74) 


(C-p ) = 2(;^)(tan 3^ + tan 3_) cos 3 

^ R ^R ^ 


89 


m 


(3.75) 


Dunham and Game suggest that the method of 
reference 88 would be applicable to a wider range of turbines 

I I T* 77 < ^ 

if the following secondary and tip cYehf^ce iQS,s correlations 
are used 


%5r« ^ 


A jasjs 



52 


X 




X 


c 


R 


0,0334 ( 

0.0334 ( 


h. 


'I 


h. 


'1 

■R 


)(■ 


cos a 


cos a 


cos 3_ 


(3.76) 

(3.77) 


Thus the sum of the secondary and tip clearance loss coeffici- 
ents of the rotor and the nozzle are given by: 




E 


2 

cos B 


c ^ 0.78 

= ^ ’■o > < fs 7c” 7 > ^ — T~ 

°R ^R ° '^®R^ R^ cos^ B 


i) 

m 


(3.78) 


cos^ a 

^c + { 7s TF T> ^ 

m 


(3.79) 


where k is the clearance and B is a constant depending upon the 
type of clearance. The value of B is 0.50 for a radial tip 
clearance and 0.25 for a shrouded blade with side clearance 
as shovm in Figures 3.3(c) and 3.3(d). 


Step (iii); The total loss coefficient for nozzle and rotors 
can be calculated as: 


Y.. = (Y ) + (Y„ + ) 

R P TST s k^ 


'R 


R 

) 

R 




R 


R 


(3.80) 


(3.81) 
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Step (iv); From the known Yjj and Yp, the iterated (improved) 




loss parameters and are obtained as; 


Y T' 
N 2 


^*1 


T 


o1 


(3.82) 


and 


E. 

1 


Y T" 
E 3 


T' 


o3 rel 


(3.83) 


where = '^3 + 


(3.84) 


It can be observed that X^^ is greater than X^^ by 
virtue of tip leakage loss in the rotor blades. Y/ith the new 
values of X„ , an improved value of the isentropic 

efficiency is calculated by making use of relation (3.7) as: 


1 


sx 


2 2 
1 p 2 1 p 


(3.85) 


1 + 


T ) 1 

o3^ ■' 


It is to be remembered that we started with a trial 
value of the nozzle loss coefficienb X^ and the isentropic 

stage efficiency n and now we obtained iterated (improved) 

s 

values of these quantities as X»j and Hg * It is desired that 

^^i ®i 

finally we should get X^^^ =’ and = n^. Both these requi- 
rements may not be met at a time. When convergence is not 
achieved, we start a new iteration taking X^, and n as the 
guess values ard finally we stop the iterative procedure whoa 
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converges within a specified accuracy. In the present work, 
the iteration procedure is terminated whenever { n - n . } ^ 

.0005- Finally the number of nozzle blades n-^j and rotor blades 

are calculated from the known value of spacing to diameter 
ratio as 


n- 


■R = 


lA 


( 3 . 86 ) 


and 


n - lA 


(3.87) 


iurther the Mach numbers corresponding to the velo- 


cities and 0^ (namely and M ) can be computed 


as ; 




V, 


( yR I3) 


1/2 


( 3 . 88 ) 


M, 


^2 TVr^T^T^ 


(3.89) 


The values of and M, are to be kept below the critical 

J 2 ■ 

value . 


3 . 5.5 Correction for Reynolds Number Effect 

The applicability of the correlations of profile 
losses given in section 3.5.4 is limited by the Reynolds number 
of the flow (Rg ) . The Reynolds number effects on* turbomachines 
are discussed in references 90 and 91 . The Reynolds number of 
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the flow should be in the region of 1X10^ to 3X10^ /with E 
defined in terms of the blade chord c, density p and the 
relative velocity at the outlet of the blade row. We can 
calculate the Reynolds number for the nozzle and rotor blades, 
(E^)jj and ^7 assuming a constant (known) value of 

viscosity y (y = ^2 = 

(R ) ^ _2_^ (3.90) 

® N y^ 

and 




R 


^3'^3^R 

V 3 


(3.91 ) 


The Reynolds number for a turbine stage E^ is taken as 
the arithmatic mean of R for the nozzle and rotor as 


\ = 


(Eg ) + (Eg ) 

® R ® R 


(3.92) 


If the value of E given by equation (3.92) differs much from 
c 

2X10 , an approximate correlation can be made to the isentropic 

efficiency by modifying it to (n_).D as follows: 
s s 

e 




= 1 - ( 


2jao 

R. 


5 0.2 

■) 


(1 - 




(3.93) 


All the steps followed in section 3.5 have been shown in the 
form of a flow chart in Figure 3.9. 
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3.6 Determination of Air Properties 

The computation of energy and woric of gases in the 

present chapter are based on the temperature -entropy diagram. 

This requires a knowledge of the values of specific heat of air 

at constant pressure, c , at various temperatures. These values 

P 

of Cp are taken from air tables and stored in the computer in 

the form of polynomials obtained by using least squares method. 

These polynomial expressions for c^ are given in Appendix B. 

An alternative approach for the calculation of energy and work 

is to make use of the enthalpy -entropy and enthalpy -temperature 

relations. Again these calculations require air tables to 

obtain the various properties. For this purpose, all the air 

92 

properties tabulated by Keenan and Kaye have been computerised 
in the form of polynomials using the technique of least squares 
and the results are reported in Appendix B. 

3.7 Computation of the Geometrical Properties of Airfoil 
Sections 

In order to evaluate the weight of a stage, the areas 
of cross section of nozzle and rotor blades are needed. More- 
over, some other geometrical properties of the airfoil like 
moments of inertia are required for computing the stresses in 
blades. Hence the computation of the geometrical properties 
of airfoil sections is considered in this section. 
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The blade angles and blade profiles were fixed in 

and lower 

section 3.5.3. The equations of upper profile y 
profile y^^^(x) were obtained as fourth order polynomials 
(T'igure 3.8) by making use of least squares fit. Let these 
curves be expressed in a general form as 


c^x"^ + c^x^ + c^x + C 2 X + c 


1 


Since the evaluation of the geometrical properties of the air 


K-/ .Up. 

foil requires y^ and y^ , these can be expressed 

9 ' 


y.2 _ a x^ + dgX^ + d^x^ + dgx^ 


4 - d x^ + d.x^ + d^x 
D 5 4 ^ 


+ d2X + d^ 


( 3 . 94 ) 


where d^ = g|, dg = 2c^c^, d^ - 2(c^c^ + c^), 

2 

dg = 2(c^c^ + ^2°^), dg = (2CgC^ + 2 c 2C^ + Cg), 

2 

d^ = 2(c^c^ + C2C^), dg = (2c^Cg + C 2 ) 

dp = 2 c^C2» ‘^'1 = °1 


QS) 


and, 

10 

y3 = (dgC5)x''2 + (dgOj + dgCpx''U (dyCj + dgC ^ + dgO,)^ 

+ (dgC^ + ^ ^ ^ 1 ^ 

V , 

+ dgCg . dgCpxS ^ ^ ^ 4^03 + d^Og + dgC 

+ (djO^ + + dgOj + dgCg + dyopx®' + (dgOg + 
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+ ^ ^ 5°2 "'' ‘^ 2°4 ^ 3^3 "*" ^ 4°2 

+ d^c^)x^ + (d^c^ + d^c^ + d^c^ + d^c^)x^ + (d^c^ + 


^ 2^4 + d^c^)x^ + (d^c^ + d2C^)x + d^c. 


( 3 . 96 ) 


Thus from the known values of c., c^, c_ , c. and c^ , the 

1 » 2 ’ 3 ’ 4 5 ’ 

2 

coefficients of the polynomials of y and y for upper and lower 
curves can be calculated. The geometrical properties of the 
airfoil section can be evaluated as follows: 


X 


(u) y(’^)(x) 


Area of cross-section of the airfoil = A =/ / dy dx 

( 3 . 97 ) 


i.e • , 


^(u) 

/ {y^'^^(x) - y^^^(x)} dx 

x'l) 


= ((o'"’ - + (o'"’ - of >)^ + (of’ - of’) 


2EI 

3 


+ (='"’ - of’)|- + (o'"’ - of ’)x} 


(u) 


( 3 . 98 ) 


whe 2 ?e the subscripts (u) and (1) represent the values corres- 
ponding to upper and lower curves respectively, further, 

, x'"’y'"’(x) 

X coordinate of the centroid = x = -rf ( [ x dx dy } 

^ x'l’ y'l’(x) 

x'"’ 

= /, x(y^^^(x) - y^^^(x) dx)} 

( 1 ) 

x^ ^ 


( 3 . 99 ) 
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^ (x) 

y coordinate of the centroid = y = -r[ J f y 

^ (1) (1)/ N 

X y (x) 

„(u) ^ ^ 

= hU f J ax] 

^(u) y(u)(^) 

Moment of inertia about x axis = I =f f y 

^ ^ (1) (1)/ N 

x" ' y" " (x) 

X ^3 3 

= 5 C/ { (x) - y^^^ (x)} dx ] 

y(“>{x) , 

Moment of inertia ^out y axis = I =f f x‘ 

yy ^ (1) Ui)r ) 

^ y (x) 

2 ;(u) 

= X^' {y^^^x) - y^^^(x) } dx j 

x<"> y(a)(x) 

Product of inertia = I = J J x y dx dy 


(n) 


= (x) - y^^^ (x)} dx] 


.( 1 ) 


If X 


( 1 ) 


(u) 


dx dy ] 


(3.100) 


^ dx dy 


( 3.101 ) 


dx dy 


(3.102) 


(3.103) 


= 0 , X 


will be equal to the chord of the blade, c. 
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3.8 Evaluation of the Objective function 

Once the blade angles and the chords of rotor and 
nozzle blades are decided, the cross sectional areas of rotor 
blade (Ag,) and nozzle blade (Ajj-) can be calculated by making use 
of the ai^ysis of section 3.7. The efficiency of the stage 
and the other blade parameters have been calculated in section 
3.5. Suitable breadth ratio is also assumed. Now we can 
calculate the objective function f(X) using equation (2.3). 

Since it is a mixed type of objective function, where the 
first part represents losses of the stage and second part 
represents the weight of rotor disc, rotor blades and nozzle 
blades, suitable values of and Kg are to be assumed depen- 
ding on the application of gas turbine . 

Since the procedure indicated in section 3.5 for 
computing the isentropic efficiency is iterative , it is necess- 
ary to see the nature of convergence of the procedure, for 
this the isentropic efficiency of a gas turbine stage is 
computed using different trial values of efficiency and con- 
vergence criteria. The results are shown in Table 3.2. It can 
be seen that the efficiency converged essentially to the same 
value for all trial efficiencies^ however, lesser number of 
iterations are required when the trial efficiency is assumed to 
be 0.9. The results obtained with two diffeiant values of e 
(the difference between the current trial value of efficiency and 
the calculated v^ue of efficierxjy in any iteration) show that 
more number of iterations are required with smaller (more 
stringent) values of e. 
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TABLE 3.1 

Mean Line, Thickness Distribution and Coordinates for Priroary 
and Secondary Turbine-blade Series 


SI. 

No. 

^c 

Primary turbine blade 
(reaction blading) A^Ky mean 
line coordinate for Cj^ =1 .0 
and thickness distribution 
in % of chord 

Secondary turbine blade 
(impulse) BiEiIi mean line 
and thickness distribution 



^c 

Mean line 
coordinate 
% of chord 

1 

{Thickness dis- 
Jtribution % 

{of chord 

f 

r 

\ 

^ t 

% of chord ; 

^t 

% of chord 

1 

0.0 

0 

0 

0 

0 

2 

3 

0.5 0.397 

1.25 0.836 

3.469 

0.336 

0.718 

2.583 

4 

2.5 

1 .428 

4.972 

1.253 

3 .282 

5 

5.0 

2.359 

6.918 

2.128 

4.041 

6 

10.0 

3.689 

9.007 

3.480 

5.007 

7 

15.0 

4.597 

9.827 

4.506 

5.592 

8 

20.0 

5.217 

10.000 

5.291 

6.995 

9 

25.0 

5.623 

9.899 

5.862 

8.063 

10 

30.0 

5.852 

9.613 

6.261 

9.025 

11 

35.0 

5.936 

9.106 

6.514 

9.727 

12 

40.0 

5.897 

8.594 

6.642 

10.000 

13 

45.0 

5.753 

7.913 

6 .652 

9.725 

14 

50.0 

5.516 

7.152 

6.551 

9.009 

15 

55.0 

5.200 

6.339 

6.342 

8.016 

16 

60.0 

4.814 

5.500 - 

6 .016 

6.908 

17 

65.0 

4.367 

4.661 

5.551 

5.848 

18 

70.0 

3.870 

3.848 

4.983 

5.000 

19 

75.0 

3.328 

3.087 

4.332 

4.312 

20 

80.0 

2.746 

2.406 

3.680 

3.624 

21 

85.0 

2.133 

1 .830 

2.824 

2:936 

22 

90.0 

1 .485 

1 .387 

1 .953 

2.248 

23 

95.0 

0.801 

1 .101 

0.948 

1 .560 

24 

100 

0.0 

0.0 

0.0 

0.0 
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TimLE 3.2 


Convergence Study of Isentropic Efficiency 


SI. 

! Trial 

For e = 

.0005 

For e = 

.00005 

• 

O 

lvalue of 
|eff id- 
le ncy 

Converged 
value of 
eff ic iency 

■■■{to. "of 
{iterations 
{required 

f 

Converged 
value of 
efficiency 

1 No . of 

1 iterations 

1 required 

ff 

1 

.60 

.92015686 

3 

.92015349 

4 

2 

.65 

.92015617 

3 

.92015351 

4 

3 

.70 

.92015556 

3 

.92015352 

4 

4 

.75 

.92015501 

3 

.92015353 

4 

5 

.80 

.92015453 

3 

.92015353 

4 

6 

.85 

.92015410 

3 

.92015410 

3 

7 

.90 

.92014484 

2 

.92015373 

3 

8 

.95 

.92015358 

3 

.92015338 

3 

9 

1 .00 

.92015306 

3 

.92015306 

3 

10 

1 .05 

.92015279 

3 

.92015279 

3 


Data: d = .432 m, ^ = .046, ^ = .04, ^ = -035, = -044, 

^2 = *358, = .960, = 272 m/sec , - .4X10^ , 

= 1100 K, m = 20 Kg, H = 250 rps , radial clearance 


is assumed. 
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I^Known data for stage design; , m, 11 | 


IT Tr 


• Known values of air properties sc^, E, Y , y 


Given a trial vector of design variables • 

2’ ^3’ ^a 


2r fl fR 
’ d ’ d ’ d ’ d ’ ^ 2 ’ ^ 3 ’ ^ 



Calculate gas conditions and blade height 
at station 2; Kor calculating Ig? ^2’ 

Po^/Pc> P 2 , Ag and A^j,- use equations (3.24) 

to ( 3 . 30 ) and for h 2 j equation (3.42) 

f ~ 


A 

'‘e from B 


Continued . . . 
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Calculate loss coefficients for rotor and nozzle; 
For calculating T.t, Y^, and Xr. use equations 

(3.80) to (3.83) M Tor X^ use (3.2) 


X 


. Calculate using equation (3.85) 

SI 


T — TT ] 

n = n . , 
s si 

T Cl f r\ -n ^ f ] 

^Is(ngi-ns)l ejf 


I-j£g Go To A 

Terminate f 


Figure 3.9. Flow Chart for Calculation of Sfficiency 
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CHAPTER 4 

EIHITE ELEMENT VIBRATION ANALYSIS OE ROTATING TIMOSHEMO BEAUS 

In finite element method, the domain of the continuum 
is discretized into a number of disjointed subdomains called 
elements defined by a set of points which are c^led nodes. 

The elements may not be connected at the nodes, but also along 
the interelement boundaries. Then the unknown field variables 
are approximated by a set of assumed functions, which are 
expressed in terms of the field variables at the nodes by 
suitable interpolation formulas. The characteristics of the 
continuum are predicted from the characteristics of the elements. 

The need for finite element method arose frpm the fact 
that (1 ) exact solutions for complex differential equations 
encountered in continuum approach have not been possible except 
for very few idealized cases and (2) with the advent of digital 
computers the approximate methods of solution like Ritz method, 
Galerkin technique and finite difference method, have come into 
vogue. Most well studied mathematically, among the approximate 
methods is the finite difference method which is mathematical 
discretization of differential equations. However, the 
method suffers from the inherent drawback that it cannot handle 
complicated boundaiy conditions, irregular change in geometry 
and material properties. Amongst the approximate techniques, 
the finite element method can be considered as one of the most 
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powerful and " vers afile physical discretizafion scheme presently 
available for the numerical solution of complex continuum 
problems using digital computers. Fonlinearity, ortho tropy , 
considerations of cutouts, irregular boundaries and complicated 
boundary conditions can all be incorporated in the analysis ■ : 
without much difficulty. 

In this chapter a new finite element is developed to 
find the natural frequencies and mode shapes of beams in 
bending -bending mode of vibration by taking into account taper, 
pre-tv/ist and rotation simultaneously. The coupling that 
exists between the flexural and torsional vibration is not 
considered. The taper and the angle of twist are assumed to 
vary linearly along the length of the beam. The element stiff- 
ness and mass matrices are derived and the effects of offset, 
rotation, pre-twist, depth and breadth taper ratios and rotary 
inertia and shear deformation on the natural frequencies are 
studied. Various special cases of beam vibration are obtained 
from the general equations derived. 

4.1 Displacement Model 

Figure 4.1(a) shows a doubly tapered, twisted beam 
element of length 1 with the nodes as 1 and 2. The breadth, 
depth and the twist of the element are assumed to be linearly 
varying along its length. The breadth and depth at the two 
nodal points are shown as b^ , t^ and b^j '^2 respectively. The 
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pretwist at the two nodes is denoted hy 0^ and . Figure 
4.1 (b) shows the nodal degrees of freedom of the element where 
bending deflection, bending slope, shear deflection and shear 
slope in the two planes are taken as the nodal degrees of 
freedom. Figure 4.1(c) shows the angle of twist 9 at any 
section z. The beam is assumed to rotate about the x-x axis at 
a speed of n radians per second. 

The total deflections of the element in the y and x 
directions at a distance z from node 1, iiamely, w(z) and v(z) 
are taken as 


w(z) = w^(z) + Wg(z) 
v(z) = "v-j^Cz) + ^3(2) 


(4.1) 


where w^(z) and ■^1 ^( 2) are the deflections due to bending in 

the yz and xz planes, respectively, and w (z ) and v (z^ are the 

deflections due to shear in the corresponding planes. 

The displacement models forwi^(z), Wg(z), ■v-^(z) and 

V (z) are assumed to be polynomials of third degree. They are 
s 

similar in nature except for the nodal constants . These 
expressions are given by: 


Wi^(z) 


u 


1 (2z^ - 31z^ +1^) + -i (31z' 


2z^) 


- (z^ - 21z^ + l^z) - -| (z^ - Iz^) 

1 ^ 
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Wg(z) 


^ (22^ - 31z^ - 1^) +1^ (31 z 2 _ 2z^) 

1-^ l2 


^ (z5 - 21^2 + l2,) - 1,2) 


v^,(z) 


^ (2z3 - 31,2 . i3) ^ :h0 (53 l,2 _ 2^3) 

1 ^ T ^ 


(z^ - 21z^ + l^z) --1^ (z^ - Iz^) 


(4.2) 



where , Ug j and represent the bending degrees of freedom 
and , Ug , and Ug are the shear degrees of freedom in yz 
plane, Ug, and represent the bending degrees of 

freedom and ^14’ ^15 ^16 shear degrees of freedom 

in xz plane . 


1- .2 Element Stiffness Matrix 


The total strain energy U of a beam of length. 1, due 
to bending and shear deformation including rotary inertia and 
-rotation effects is given by; 


U 


1 El 

/ [ ( 

O 3Z 


2 


2 

3 W^ 


3Z 


3^ V, 


3Z 


o + El 
2 J7 


3^v, 2 
(-^) } 
3Z 
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1 gw, 3W 2 

+ — + (iw) 1] + 5 / r(^)(-^ + t/) d^ 


1 


Sv_ 2 


1"^ ov, OV il 1 

2 / dz - / p^(z) (w + w ) d 

O n 


/ P^(z) (v^ + Vg) dz 


(4.5) 


where E = Young’s modulus, I = moment of inertia about xx aocis , 

I = moment of inertia about yy axis, I = moment of inertia 
yy ' xy 

about xy axis, p = shear coefficient, G = shear modulus and 


E(z) = /' 


1+e 


m 5 d? 


'V 


E+Z-+Z 

w 


P^A P 

[(1 + e)2 _ (E + 


+ z 


)2] 


(4.4) 


where the mass per unit length (m) is assumed to be a constant 


for simplicity. 



Pw(2) 

= ” 11 ,^ K + ™s ’ 

(4-5) 

<! 


(4.6) 


where e is the offset p is the mass density and z is the 

in 0 

distance of the first node of the element from the root of the 
beam as shown in Figure 4.1(d). 

As the cross section of the element changes with z 
and as the element is twisted, the cross sectional area A, and 
the moments of inertia I , I and I will be functions of z; 

JUi. yy JLy 
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A(z) = b(z) t(z) = + (b^ - b^) + (tg - t^) f > 

= ^ + Cglz + c^l^)' (4.7) 

where 


Cl = 

(b2 - b^ ) (t2 - 

”4 




o 

ro 

11 

bi (t2 - t^ ) + t 

^ (b2 - b^ ) 1 

> » 


(4.8) 

Oj = 







= I^,^, cos®© 

T • 2^ 7 

+ I , , sin © : 
11 \ 

1 




= ly.y. =03^6 

+ I^,^, sto®© ' 

> 

(4.9) 


= (I . t - I 
x*x’ y 

\ sin 2© 

,y,; 2 

*> 




where x'x' and j'j' are the axes inclined at an angle ©, the 
angle of twist, at any point in the element, to the original 
axes XX and yy as shown in Figure 4.1(c). The value of I 

X J 

0 ard the values of I„,vi and I can be computed as; 

XX 


T 

“X 'x 


(z) 


b(z) t^(z) 
12 


where 


121 


^{aiZ^ 


T 3 , n2 2 
+ a^lz + a^l z 
2 i 

4 

+ a^l z + a^l } 


( 4 . 10 ) 


a^ = (bg - h^ ) (tg - )^ 

^2 = ^1^'^2 " ^1^^ + 3(b2 - ■b-,)(t2 - -t^)^ t^ 
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= 3 _ b^)(t2 - t^)t2 


a^ = 3lD^t^(t2 - t^) + (b^ - b^)t5 

®5 = ” 1^1 


(4.11) 


T ( 7 ) - • 'fa(z) .b^(z) 

7'7'^ - 12 


^ {d^z'^ + dglz^ + 


121 


+ d^l^z + d^l"^} 


( 4 . 12 ) 


where 


di = (t^ - t^)(b2 - b^)^ 

d2 = t^(b2 - b^)5 + 3(t2 - t^)(b2 - t-,)^b^ 

d^ = 5 {t^b^ (b2 - b^)2 + (t^ - t^)(b2 - b^ )b^ } (4.13) 

d^ = 3t^b^(b2 - b^) + (^2 - 't'l 

45 = t^b3 

By substituting the expressions of w^, w^, v^, , A, 

I , I and I from equations (4.2), (4.7) and (4.9) in 
XX ’ xy 77 

equation (4.3)? the strain energy U can be expressed as: 


U = i u ’^[K] 3 (4.14) 

where u is the vector of nodal displacements u^ , U2 ,...,u^g, 
and [kJ is the elemental stiffness matrix of order 16. 

Denoting the integrals 
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1 

r 

CM 

CO 

2 

/ 


•) dz = [ u 

0 

3Z 


1 

El (t-k^ 

2 

/ 

o 

■) dz = [ u, 

1 

3W 2 


/ 

0 

VAG(— ^) 
9Z 

dz = [u^ ^ 

a 

2 

9 W, 

2 

/ 

0 

El ( — 5^ 

) (— 2^)dz 
9Z 

1 

3W 

2 

J 


dz = [ u^ 

0 


and 


1 

I 

2 P^A (2^(w- 

^)dz = [ u^ 

0 


T 


T 


,T 


(4.16) 


(4.17) 


iT 


T 


(4.18) 


(4.19) 


, '^4^^ [5’K][u^ Ug ] (4.20) 

The element stiffness matrix can be expressed as 
[AK]+[EK]-[T’K] [EK]-[EK:] [DK] [O] 

[Eig -[EK] [ ciq +[ek] -[ek] [o] [ o] 

[DK] [ 0 ] [bk] +[ek] -[eiO [ ekO -[ek] 

[o] [o] [ek]~[ek] [ck]+Cek] -[ek] 

( 4 . 21 )' 

where [AK], [BK] , [OK], [ DK] , [ EK] and [EK] are symmetric 
matrices of order 4 and their elements are formulated in 
Appendix 0. [0] is a null matrix of order 4. 


[k] = 
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4.3 Element Mass Matrix 

The kinetic energy of the element T including the 
effects of shear deformation and rotary inertia is gi'ven by 


T = J ^ [A(— ^ + -^) + A(-4 + -^) + I ( ^ 

9t at az at 


o 


a t 


) 


a t 


^2 
a W-. 


a^v. 


+ 21 ( 2_)(- 

^ az at az at 


2 

aw, 2 

-) + I ( — ] az 

az at 


( 4 - 22 ) 


By defining 


1 aw, 2 

/ ( z) = 


m 


at 


1 9^w, 2 

/ p„,I ( dz 

J iric ^ 


m XX ^ ^ . 

0 3 Z 9 t 


1 aV 2 

f P T ( ) dz 

o “S^s^azat 


[ 1^2 U3 U^] ^ [M4][u^ Ug S V ’ 

= [u,, 1^2 % ^2 S ’^4^’ 

(4.24) 

= [ Ug ^11 2 ^ ^ ^ ^9 ^10 ^11 ^ 12 ^ 

(4.25) 


ard 

1 a^w 

f = [*1 *2 a , ^0 *11 * 12 ^ 

0 •' az at az at 

(4.26) 

where Uj_ denotes the time derivative of the nodal displacement 
u^, i = 1, 2,... ,16, the kinetic energy of the element can be 


expressed as 
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T = ^ u ’’ [H] u 


(4,27) 


where M is the mass matrix given by 


Ck] 

16 X 16 


[i\lO +[B^ 

[AM] 

[DM] 

[0] 

[. M] 

[M] 

[am] 

[ 0] 

[ m] 

[m] 

[M +[ CI^ 

[am] 

[0] 

[0] 

[ Aid 

[am] 


(4.28) 


and £ AM] , [ BM] , [ CM ] and [ DM] are symmetric matrices of order 
4 whose elements are defined in Appendix C . 


Boundary Conditions: 

The following boundary conditions are to be applied 


depending on 

the type of 

end 

conditions; 




3W 


3Y 




Dree end : 

= 0 

9z 

and 

S 

dZ 

= 0 


(4.29) 

Clamped end : 

w„ = 0, w, 
s ' b 

= 0, 

\ = o> 

li 

O 

9w 

— ^ _ 0 

9 Z 




= 0 




(4.30) 

Hinged end; 

Wg = 0, 

= 0, 

V = 0 
s 

and v^ 

= 0 

(4.31) 


4.4 Special Cases 

The various special cases of the beam vibration 
problem can be solved by applying one or more of the following 


four conditions; 
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(a) For non-rotating beams? 

= 0 which results in [ EK] = [ EK] = [ 0 ] (4.32) 

(b) For uniform beams; 

By setting = b^ and t^ = t^ , one obtains 


= °2 

= 0 

and c„ = 
3 

hb 



Q, ^ ~ 

— Siry 

= a^ = 0 

and 

^5 = 

(4.33) 

CM 

!1 

11 

11 

11 

O 

and 

b = b'=1 



(c) For neglecting the effect of shear deformation; 

w = V = 0 so that equations (l) and (2) become 
s s 

w(z) = w^(z) and (4.34) 

Due to this, the order of [Kj and [m] matrices reduces 
from 16 to 8, 


(d) For beams without pre-twist; 

In this case there will be no coupling between the moment 
of inertia terms and one obtains; 


I = I , , 

XX x'x' 

I = I , , 

77 7'7 

I =0 

xy 

For vibration in yz plane , ~ 

For vibration in xz plane, ^ 


(4.35) 




Pew of the special cases of practical interest are 
discussed below. 


4.4.1 Classical Tapered, Twisted and Rotating Beam 

This case is obtained by neglecting shear deformation 
in the general case. When the conditions of equation (4.34) 
are applied, the [K] ani [M] matrices reduce to; 

, , r [M] + [EK] [PK] [m] 

[KJ = : ^ 

[ DKj [BK] + [ EK] - [PK] 

and 




[AM] + [BM] 
[ DM] 


[ DM] 

£ Aivi ] + [ cm] 


(4.37) 


4.4.2 Mon-rotating , Tapered and Twisted Beam with Shear 
Deformation and Rotary Inertia 

This case can be obtained by applying the condition 
of equation (4.32) to the most general case. Here the mass 
matrioc [M] remains same as one given in equation (4.28) while 
the stiffness matrix [K] reduces to the following 

[AK] [0] [DK] [0] 

[0] [CK] [0] [ 0 ] 

[DK] [0 ] [BK] [ 0 ] 

[0] [0] [0] fCK] 
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4.4 *3 Classical Tapered and Twisted Non— rotating Beam 


By applying the non-rotating condition of equation 
(4.32) along with equation (4.34), the stiffness and mass 
matrices reduce to order 8 and are given by 


and 


Ck] 


[AK] [IK] 

[DK] [BK] 


[M] 


[m] + [bm] [dm] 

[dm] [m] + [cm] 


(4.39) 


(4.4€) 


4.4.4 Non-rotating, Untwisted and Tapered Beam with Shear 
Deformation and Rotary Inertia 


In this case the non-rotating condition of equation 
(4.32), the uniform beam condition of equation (4.33) and the 
no-twist condition of equation (4.35) are applied simultane- 
ously, The beam may vibrate either in yz plane or in xz plane. 

In the special case of a non-rotating, untwisted beam 
vibrating in the yz plane the element stiffness and mass 
matrices will be of order 8 and can be expressed as follows t 


and 


“[A] [0]^: 

:^[o] [B]_ 

’’[c] + [D] [C] 

^ [C]’ [C] 


(4.41 ) 


[m] ■ = 


(4.42) 
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where [A], [B], [C] and [D] are symmetric matrices of order 4- 
whose elements are given in Appendix D. 

4.4.5 Classical Fon-rotating, Untwisted and Tapered Beam 

If shear deformation is neglected from the case 
discussed in section 4.4.4, the following simplified fom of 
stiffness and mass matrices of classical, non-rotating, 
untwisted and tapered beam will be obtained; 


i — 1 

= [A] 

(4.43) 

[M] 

= [C ] + [D] 

(4.44) 


where [A], [C] and [ D] are symmetric matrices of order 4 whose 
elements are given in Appendix D. 

4.4.6 Non-rotating, Untwisted and Uniform Beam with Shear 
Deformation and Rotary Inertia 

If the conditions (a), (b) and (d) are applied one 
gets the following expressions for [K] and [M] for non-rotating 
uniform beams without pre-twist but with a consideration of 
rotary inertia and shear deformation effects (for vibrations in 
y-z plane ) ; 
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El 


K = 


XX 


1 


-12 ,61 

-61 

0 

0 

0 

0 

12 61 

61 

0 

0 

0 

0 

CM 

H 

21^ 

0 

0 

0 

0 


4U 

0 

0 

0 

0 



36J 

-36J 

-31J 

-31J 




36J 

31J 

31J 





p 

41 J 

-l^J 

mmetric 





C\j 

H 


(4.45) 


M 



156+36P 54-36P 

221-31P 

131-31P 

156 

1 

54 

221 

131 : 


156+36P 

-131+31P 

221+3 IP 

54 

156 

-131 

221 



2 2— 
41 +41 P 

2 2— 
-31 -IP 

221 

-131 

41^ 

-31^ ■ 




2 2— 
41 -41 P 

131 

221 

-31^ 

41^' 

420 ' 




156 

54 

221 

131 






156 

-131 

221 ; 


Symmetric 





41^ 

-31^ 


Ki- , 






,41^^ 


where J = 


30ll 


(4.46) 


XX 


and 


UI 


P = 


XX 




4.4.7 Classical, Fon-rotating, Untwisted and Uniform Beam 

Equations (4.45) and (4.46) further reduce to the 
following well-known equations if the effects of shear deforma- 
tion and rotary inertia are neglected. 
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[K] 


[M ] = 



12 -12 

-61 

-61 

El 

XX 

12 

61 

61 

P 


41^ 

CM 

H 

CM 


Symmetric 


41^ 

_ 


1 1 56 54 

-221 

131 


156 

-131 

221 

420 


41^ 

-31^ 


Symmetric 


41^ 


(4.47) 


(4.48) 


4.5 Numerical Results: 

The element stiffness and mass matrices developed are 
used for the dynamic analysis of cantilever beams. By using 
the standard procedures of structural analysis, the eigen value 
problem can be stated as: 

([K] - 0)^ [M]) ^ = 0 (4.49) 

where [K] and [Mj denote the stiffness and mass matrices of 
the structure respectively, U indicates nodal displacement 
vector of the structure, and w is the natural frequency of 
vibration. 


4.5.1 Convergence Study 

A study of the convergence properties of the element 
is made by taking special case of a unifoim beam discussed by 
Thomas ard Abbas^®. The results are shown in Table 4.1,. and 
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the convergence can he seen to be quite good. It can be 
observed that the eigen values converge to the exact values 
even v/ith four elements. The convergence of the natural frequ- 
encies of a rotating pre twisted doubly tapered cantilever beam 
has also been studied by considering the effects of rotary 
inertia and shear deformation and the results are shown in 
Table 4 •2. In this case also it is found that reasonably 
accurate results can be obtained even by using four finite 
elements . 

4.5.2 Rotating Beam Results 

The effects of shear deformation and depth taper 
ratio on the natural frequencies of a rotating twisted beam 
are shown in Figure 4.2 for a beam of length 0.254 m, offset 
zero, depth at root 0.00865 m, breadth at root 0.0173 m, twist 
45°, rotation 100 r.p.m. and breadth taper ratio 3. The 
material properties of the beam are taken same as those given 
in Table 4.2. The effect of shear deformation is found to 
reduce the frequencies at higher modes while at lower nodes the 
results are nearly unaffected. There is an increase in the 
frequencies of vibration with an increase in the depth taper 
ratio in the first, second and fourth modes while a decrease 
has been observed in the case of third mode (vibration in a 
perpeMicular plane). Figure 4.3 shows the variation of 
natural frequencies with breadth taper ratio. In this case the 
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data is same as in the case of Figure 4.2. 

In Figures 4.4 and 4.5 the variation of frequency 
ratio with rotation and pre-twist is studied. It can he seen 
that frequency ratio changes slightly with the rotation but 
appreciably with the twist. At higher modes (in Figure 4.5) 
the effect of twist can be seen to be more pronounced. It is 
also observed that the frequency ratio increases with an 
increase in the twist in the case of first and third modes while 
it decreases with an increase of the twist in the case of 
second and fourth modes of vibration. 

In Figure 4.6 the effect of offset is studied for a 
twisted blade having 60° twist with other data same as that of 
Figure 4.4. The frequency values are shown in Table 4.5. It 
is observed that an increase in offset changes the frequency 
ratio more at higher values of rotation. The frequency ratio 
has been found to increase with an increase in the offset . 

4.5*3 Non-rotating and Twisted Beam Results 

Figure 4.7 shows a comparison of the results given 
by the finite element method with those reported by Rosard for 
a uniform beam of 0.0254 m x 0.00635 m cross section and 
0.2794 m length for 0° to 30° twist. It can be seen that 

results are quite comparable. 

Figures 4.8 and 4.9 represent the frequency ratio 
of an uniform beam and a twisted beam for first four modes 



91 


when the angle of tv^ist is varied from 0° to 90°. The length 
of the beam is taken as 0.254 m and the cross section as 
0.076 X 0.038 times the length of beam. The results a.xe shov/n 
for Timoshenko beams where shear deformation effects are 
considered and for beams where shear deformation effects are 
neglected. The frequency values are also shown in Tables 4.4 
and 4.5* It can be seen that the shear deformation effect is 
comparatively more at higher modes of vibration and also 
frequency ratio changes at a higher rate with an increase in 
the twist angle at higher modes. The present results are 

p p 

found to agree well with those given by Mabie and Rogers 
(Table 4 .4 ) . 

Figures 4.10 and 4.11 give the variation of modal 
frequencies with breadth taper ratio for beams having 0°, 30°, 
60° and 90° twist with constant depth taper ratio while 
Figures 4.12 and 4.13 show similar variations for beam with 
constant breadth taper ratio and varying depth taper ratios. 
The beam dimensions are same as used for Figures 4.8 and 4.9. 
Again the effects of breadth and depth tapers are seen to be 
pronounced at higher modes of vibration. Here also the effect 
of shear deformation is seen to reduce the modal frequencies 
at higher rate at higher modes of vibration in all the cases. 

The effect of cross sectional dimensions for the same 
area of cross section on the natural frequencies is also 
studied and the results are shown in Table 4.6. The present 
results are found to agree well with those available in the 
literature . * 
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4.5.4 Non-rotating and Untwisted Beam Results 

Tables 4.7 and 4.8 show a comparison of the natural 
frequencies of a non-rotating untv/isted tapered beam for 
various combinations of depth and breadth taper ratios. Sioc 
finite elements are used to model the beam. It is observed 
that for constant depth taper ratio the frequency ratio of all 
the four modes increases with breadth taper ratio while for 
constant breadth taper ratio the frequency ratio decreases 
for the first mode and increases for the second, third and 
fourth modes with an increase in the depth taper ratios. The 
shear deformation effects reduce the frequency of modal 
vibration. The present results can be seen to compare well 
with those reported by Mabie and Rogers. 

4.6 Conclusion 

The finite element procedure developed for the eigen 
value analysis of rotating, doubly tapered and twisted 
Timoshenko beams has been found to give reasonably accurate 
results even with four finite elements. The effects of breadth 
and depth taper ratios, twist angle, shear deformation, offset 
and rotation on the natural frequencies of vibration of 
cantilever beams are found. The element developed is expected 
to be useful for the dynamic analysis of blades of roto- 
dynamic machines. 
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TABLE 4.1 

Natural Frequencies (Hz) of an Untwisted Uniform Beam 


Number of 
elements 

First mode 

Second mode 

Third mode 

Fourth mode 

1 

849.3 

5005.0 

12278.4 

23240.1 

2 

846.1 

4012.2 

9531 .6 

15877.0 

3 

845.9 

3996.3 

8895.3 

14410.0 

4 

845.8 

3991 .8 

8861 .4 

13909.6 

5 

845.8 

3990.4 

8846.6 

13865.5 

6 

845.8 

3989.8 

8840.7 

13843.1 

7 

845 .8 

3989.7 

8838.1 

13832.4 

8 

845.8 

3989.5 

8836.8 

13827.1 

Exact 

845.8 

3988.9 

8834.2 

13818.1 


Data: Length of beam = .254 breadth - .0762 m, 

depth = .08 fl2 X length of beam, E = 2.07 x 10 

* 2 ^ 

N/m^, Gr = 3/8 E, mass density = 800 kg/m" , 
shear coefficient = 2/3. 
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TABLE 4.2 

Convergence Study of Natural Frequencies (Hz) of a 
Tapered, Twisted and Rotating Beam 


Number of 
elements 

First mode 

Second mode 

Third mode 

Fourth mode 

1 

304.30 

1 191 .66 

2327.07 

4552.23 

2 

296.22 

1161.70 

1779.50 ' 

4106.83 

3 

295.03 

1155.97 

1746.55 

3747.04 

4 

294.85 

1154.94 

1741 .39 

3697.05 

5 

294.78 

1154.67 

1739.99 

3689.82 

6 

294.78 

1154.58 

1739.40 

3683.98 

7 

294.78 

1154.54 

1739.25 

3683.91 

8 

294.78 

1154.50 

1739.10 

3683.85 


Data: Length of beam = 0.1524 m, breadth at root = .0254 m, 
depth at -root = .0046 m, depth taper ratio = 2.29, 
breadth taper ratio = 2.56, twist angle = 45 , 
shear coefficient = .853, E = 2.07 x 10 N/m , 

3 

G = E/2.6, offset =0, mass density = 800 hg/m , 
rotational speed = 250 revolutions/second. 
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TABLE 4.3 

'E^feci; of Offset on Natural Erequencies (Hz) of a Rotating Twisted 
Beam 


Mode 

shape 

; Twist 
» 


Offset 

= 0 

f 

f 

J 

f 

Offset = 

0.0254 

in 

|rps = o; 

« * 

100 i 

• 

200 ; 

» 

300 i 

- - f 

0 : 

f 

100 ; 

1 

200 j 

300 


0° 

286.4 

271 .9 

222 .6 

91 .0 

286.4 

276.7 

244.8 

177.5 

T 

30° 

299.5 

285.3 

237.5 

119.7 

299.5 

289.8 

258.3 

193.5 


60° 

323.2 

309.7 

264.8 

164.3 

323.2 

313.8 

283.5 

223.4 


90° 

347.1 

334.3 

292.4 

203.9 

347.1 

338.1 

309.2 

253.6 


0° 

901 .3 

897.0 

883.0 

861 .0 

901 .0 

898.0 

889.0 

874.0 

TT 

30° 

890.0 

886.0 

872.0 

850.0 

890.0 

887.0 

878.0 

863.0 

X X 

60° 

861 .4 

857.0 

844.0 

821 .0 

861 .4 

859.0 

850.0 

835.0 


90° 

810.6 

806.0 

793.0 

770.0 

810.6 

808.0 

800.0 

786.0 


0° 

1789.0 

1802.0 

1838.0 

1898.0 

1790.0 

1806.0 

1854.0 

1931 .0 

TT T 

30° 

2081 .0 

2092.0 

2124.0 

2176.0 

2081 .0 

2096.0 

2138.0 

2206.0 

111 

60° 

2569.0 

2578.0 

2604.0 

2647.0 

2569.0 

2581 .0 

2615.0 

2672.0 


90° 

2993.0 

3000 .0 

3022.0 

3059.0 

2993.0 

3002.0 

3032.0 

3080.0 


0° 

5478.0 

5482.0 

5493.0 

5512.0 

5478.0 

5483.0 

5498.0 

5523.0 

Ttr 

30° 

5219.0 

5224.0 

5241 .0 

5267.0 

5219.0 

5226.0 

5247.0 

5282.0 

IV 

60° 

4799.0 

4805.0 

4820.0 

4846.0 

4799.0 

4806 .0 

4826.0 

4860.0 


90° 

4|39.0 

4045.0 

4064.0 

4094.0 

4039.0 

4047.0 

4071 .0 

4111.0 


Continued. . . 
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Table 4.3 (Continued) 


Mode 


Offset = 

t - 

0.0508 

m 

! Offset = 

7 

0.0762 

m 

shape 

J Twist 

'EPS = 0 
! 

100 

! 

200 

1 

i 300 

t 

! 0 
t 

100 

1 

200 

i 300 


0° 

286.4 

281 .4 

265.2 

233.9 

286.4 

286.0 

284.0 

279.0 

T 

30° 

299.5 

294.0 

277.5 

246.0 

299.5 

298.6 

295.5 

289.1 

X 

60° 

323.2 

317.8 

300.9 

269.8 

323.2 

321 .8 

317.4 

309.3 


90° 

347.1 

341 .8 

325.2 

295.1 

347.1 

345.5 

340.4 

331 .4 


0° 

901.3. 

900.0 

895.0 

887.0 

901.0 

901 .0 

901.0 

900.0 

II 

30° 

890 0 

889.0 

884.0 

877.0 

890.0 

890.0 

890.0 

890.0 

60° 

861.4 

860.0 

856.0 

850.0 

861 .0 

862.0 

863.0 

864.0 


90° 

810.6 

810.0 

807.0 

801 .4 

810.6 

811 .0 

813.0 

817.0 


0° 

1789.0 

1810.0 

1869.0 

1964.0 

1789.0 

1814.0 

1884.0 

1997.0 

III 

30° 

2081 .0 

2099.0 

2152.0 

2236.0 

2081 .0 

2103.0 

2165.0 

2266.0 

60° 

2569.0 

2583.0 

2627.0 

2698.0 

2569.0 

2586.0 

2638.0 

2723.0 


90° 

2993.0 

3005.0 

3042.0 

3102.0 

2993.0 

3007.0 

3052.0 

3123.0 


0° 

5478.0 

5484.0 

5503.0 

5533.0 

5478.0 

5485.0 

5507.0 

5544.0 

IV 
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Data; Length of beam = 0.1524 m, breadth at root = .0254 m, 
depth at root = 0.00103 m, depth taper ratio = 1 , 
breadth taper ratio = 1, shear coefficient = .833? 
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E = 2.07 X 10^^ N/m^, G = E/2.6, mass density = 800 kg/m , 
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MLE 4.6 

Effsci; of fixed ind Cross Sectional Dimensions on Daiural 
frequencies (Hz) of lapered Twisted Beam with Shear 
Deformation Effects having Constant Area of Cross Section 
at fixed End 


Ratio of breadth 
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to jPre- ; 


Mode 

shapes 


thickness at fized ;twistl 
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4786 .0 

2.0/1 .5 

twist) 
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5982.5 


2 

Data: Length of beam = .254 m. Constant root area = .000762 m , 
Depth taper ratio = 3, Breadth taper ratio = 2, 

Shear coefficient = .853 > Mass density = 800 Kg/nP 
E = 2.07 X lo'^^ N/m^, G = 3/8 E. 
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CHAPQIER 5 

EVALUATION OP CONSTRAINTS 


The optimization problem has been formulated in 
Chapter 2 and the method of calculating the objective function 
was discussed in Chapter 3. In Chapter 4, the vibration 
analysis of tapered pretwisted rectangular beams has been 
considered using the finite element method. The method of 
evaluating the constraints of the optimization problem is 
discussed in this chapter so that the numerical techniques of 
optimization discussed in Chapter 6 can be applied for solving 
the problem. 

5.1 Blade Idealization 

Generally the turbomachine blades are idealized as 
cantilever beams of rectangular cross section. Refinements in 
idealization can be made by including the effects of taper, 
pretwist, rotation, shear deformation and rotaiy inertia as 
described in Chapter 4, In order to apply the analysis 
procedure of Chapton 4j the airfoil section has to be replaced 
by an equivalent rectangular section with pretwist . A method 
of replacing a rotor blade of airfoil section by a doubly 
tapered and pretwisted Cantilever beam of rectangular cross 
section is described in the following subsections. 



1 ( 


5.1.1 Equivalent Rectangular Section of an Airfoil Section 

Prom the known values of area A and moments of 
inertia ^yy (computed in section 3.7)> the 

moments of inertia I , I and I with respect to x and y 

XX xy yy 

axes passing through the centroid of the airfoil section can 
be computed using the parallel axes relations as; 


I__ = 

ioc 

I 

XX 

C\J 

1 

(5.1) 

II 

H 

I 

yy 

- A5E2 

(5.2) 

I 

xy 

I 

xy 

- A 2^ 

(5.3) 


We then define a rectangular section having tv/ist ©, 
breadth b^^ and height t^^ to be equivalent to the given air- 
foil section by equating their areas, and moments of inertia 
about X and y axes (Pigure 5.1)' Thus 

Area of the rectangle = A = b^^ t^^ (5.4) 

Moment of inertia of the rectangle about x axis = 


XX 


12 ^eq ^eq 


^ cos^& + -t.,. b:'.. sin"© 


12 eq eq 


(5.5) 


Moment of inertia of rectangle about y axis - 


1 


77 


12 ^eq ^eq 


5 cos^e + h b„. sin^© 


12 eq eq 


(5.6) 
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The three unknown parameters ©, b and t can be 

eq eq 

found b 7 solving the three equations (5.4) to (5-6) simultane- 
ously. Here it is assumed that I is negligible compared to 

icy 

I and I . Prom equations (5.4) and (5.5), one gets 

XX yy 


121 


XX 


poop 9 9 9 9 9 

— = t „ COS © b ^ sin 9 = t~„ cos 9 + b - b „ cos 9 

A eq eq eq eq eq 

(5.7) 


i.e . 


2 

cos 9 


12 I 


XX 


b^ A 
eq 


A(t 


eq 


b^ ) 

eq^ 


(5.8) 


Similarly from equations (5*4) and (5.6), 

2 


sin 9 


2„ yy 1 


A(t 


eq 




(5.9) 


Equations (5.8) and (5.9) give 


12 I__ - A b 
XX 


eq 


eq 


A(t 


’eq 


b2 ) 

eq^ 


12 I - A b 

yjj 

^^’'"eq ” ^eq^ 


= 1 


( 5 . 10 ) 


By substituting the value of b^^ from equation (5.4), equation 
(5.10) gives 


.4 _ 12 

'eq ~ A 


(1. + I—) ^la + ^ 

XX yy 


0 


(5.11 ) 


where roots are given by 
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XX 


+ _) i 

jy 


A xx yy 


)^ - ( 5 . 12 ) 


and 


t) 


eq 


A 

^eq 


(5.13) 


Finally the twist can he obtained from equations (5.8) and 
(5.9) as 


, ^ sinG 

tan© — cos© 


( 021 - 

yy 


0.5 


(5.14) 


Notice that equations (5.12) and (5.13) gi'^Q tv/o 
equlTOlent rectangular sections whose orientations are perpen- 
dicular to each other. Out of these, we select the one hartng 
a larger value of h^^^ compared to t^^ as the desired equivalent 

section. 

T Tt- c; 1 -7 axis will be in the radially out- 

In Figure 5.i> z-axis 

ward direction, r-accis will he along the chord direction and 

rya ares form a right handed coordinate system. The origin of 

4 - n-f +He cxoss S0cl3ioii of fiio 
xyz system will "be at the cen r 

hlade at the root. Since the airfoil section changes from 
point to point along the length of the hlade (s-aris), e 
value Of . Which defines the directions of degrees of freedom 
With respect to r'-aris will also change from point to poxn . 
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5.1.2 Equivalent Doubly Tapered and Twisted Rectangular 
Cantilever Beam of a Rotor Blade 


Generally turbine blades having airfoil shapes are 
tapered from root to tip as shown in Figure 5.2. The taper is 
such that the value of (area of tip/area of root) lies betv^een 
1/4 and l/3. In this work the taper in chord or breadth (g) 
has been fixed as 2 and the taper in thickness (a) has been 
assumed as 1.5. The chord leugtos at root and tip are calculated 
in terms of the known chord length at mean radius (Cjj) as; 



root 


2b 

(1 + ?T 


(5.15) 


(Cp) 

^ tip 


(cj^) /e 

" root 


( 5 . 16 ) 


The hlade angles Sgroot’ ^3root 

fc.. at tip are eaready calculated from the free vortex 
'■^3 tip 

analysis. 

Once the values of (cg)root’ ^2root 3root 

known, we can get the equivalent rectangular section of 

.V. ^ 1 and twist 6^^ . at root by 

dimensions (hQq)poot’ ^^eq'^root ^ 

R 1 1 . Similarly using the 
using the procedure of section 5. 

R we can obtain an equivalent 

values (Cj^ ^tip’ ®2tip ^3tip 

rectangulax section of diniensions tip^ 

twlBt S.. at tip. Hence the pretwist of equivalent blade 

tip 


Will be 


0 ®4--f TV 




'root 


(5.17) 
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Thus the actual rotor blade is idealized as a 
doubly tapered cantilever beam of length hj^ and pretwist 0 
whose end dimensions are and at the root 

and (b )j.- and (t . at the tip. If the airfoil blade 

ciixv^ tip eq tip 

itself has an original pretwist of the total pretv'/ist of 
the blade v/ould be (© + ) . 

5,2 Hot or Blade Stresses 

After obtaining the idealized (equivalent) blade 
section we have to verify whether the stage design is consistent 
with the permissible stresses in the rotor blades. Bor this we 
assume that simple approximate methods are adequate for pre- 
liminary design calculations. The final design has to be checked 
by laying out the blade cross section at several radii between 
root and tip, and performing an accurate stress analysis on the 
lines indicated by Sternlicht^^ . Mong the various types of 
stresses induced in rotor blades, only centrifugal stress, gas 
bending stress and stress due to pressure force are important 
and hence are considered in this work. 

5.2.1 Centrifugal Tensile Stress 

When the rotational speed ot the rotor (H revolutions/ 
second) is specified, the allowable centrifugal tensile stress, 
places a limit on the annulus area but It does not 
affect the chord length, of the blade. The maximum value of 
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the centrifugal stress occurs at the root of the blade and is 
given by 


^ '^ct^ 


m 


2 r 

tip 


max 


A 


'root r 


A^r dr 


'rooT: 


(5.18) 


where is the mass density of the blade material, .q is the 
angular velocity of the blade (rad/sec), A^ is the cross 
sectional area of the blade at any radius r and is the 

cross sectional area of the blade at root. In practice the 
integration of equation (5.18) is performed either graphically 
or numerically. However, if the blade is of uniform cross 
section, equation (5.18) reduces to 




max 


’’..a l , (3,2 - ) = 

2 ^^tip root 


2 

p „ fi A _ 
m gn 

2 TT 


2TrN^ '^an 


(5.19) 


where ^root radii of tip and root of the 

blade respectively. A „ is the annulus area of the turbine at 
mean height of rotor blade. A rotor blade is usually tapered 
in chord and thickness. For preliminary design calculations 
it is sufficiently accurate (and on the safe side) to assume 
that the assumed taper reduces the stress to 2/5 of the value 
for an untapered blade . Thus the approximate value of the 
maximum centrifugal stress induced in the rotor blade can be 
taken as 
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(a 


ct' 


max 


i 

3 


p A 
' m an 


( 5 . 20 ) 


Note that in a multistage turbine, the maximum 
centrifugal stress will occur at the root of the blade of the 
last stage where the fluid density is least and the aruiulus 
area, therefore, is largcsst. 

The centrifugal stresses can also be calculated 
using the finite element technique. For this, two additional 
degrees of freedom (in axial direction) are required in each 
element, which will increase the size of the elemental stiffness 
and mass matrices . The total number of degrees of freedoms 
for the blade will also increase which means some more compu- 
tational time. Hence the simpler procedure stated above has 
been employed in calculating the centrifugal stresses. 


5.2.2 Stresses Due to las Bending and Pressure Force 

The force arising from the change in the angular 
momentum of the gas in tangential direction, which produces the 
useful torque, also produces a gas bending moment about the 
axial direction, namely, as shoY/n in Figince 5 . 5 - ihere 
may be a change of momentum in axial direction (when 0^^ ^^^2^ 

and with reaction blading there will be a pressure force 

\ ( 5 . 21 ) 


(Po “■ Px^ 




per blade of the rotor where d is the diameter of rotor, h^^ is 
the height of rotor, np is the number of rotor blades and P2 
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sUd. QJC& ijiiG piBSSui'ss. Tiiis prBSsui'B fopce will givs piss 
to a gas bending moment about the tangential direction-. By 
resolving these bending moments into components parallel to the 
principal axes (X - Y) of the blade cross section, the maximum 
stresses can be calculated. The moments along X and Y are 
given by 


Along X 

Along Y 
so that 


M = ll cos(© 
X 


M_ = M cos(© 

y 3 . 


a ) + M sin(© - a ) 

o d. o 

a ) - M sin(© - a ) 
s w , s 


stress induced 


Y 

XX ^ 



( 5 . 22 ) 


Equation (5.22) assumes that M_ and M_ are constant along the 

X Y 

length of the blade. The stress induced in the blade can be 
computed more accurately if the variation of M_ and along 
the length of the blade is considered. This can be done without 
much difficulty using the finite element formulation. 

For this, a twisted and tapered blade has to be 
divided into strips of height 5 h and the bending moments 
calculated from the average force acting on each strip. For 
simplicity, these strips can be taken to be same as the finite 
elements used for the blade analysis. In the present work, 
the stresses due to gas bending and pressure forces have been 
calculated by using the finite element technique discussed in 



125 


Cha.p'ttii’ 4 3-S iij will nob I'ef^.uipo any addilional dognoes of 
fnaodoiQ. pon ©lonionf. Pox' calc tils, "ting "bh© sinsssssj foxcos ax6 
applied at the nodal points in the directions of the nodal 
degrees of freedom and then the deflections due to loading are 
obtained by applying standard techniq.ues of matrix structural 
analysis. Once the deflections are available, the stresses at 
any section can be calculated by using the relations from the 
theory of structures. The procedure adopted in finding the 
vector of nodal loads is given below. 


(i) Forces due to gas bending: 


The quantity (C^ + 0^ ) represents the change in 

2 3 

the whirl (or tangential) component of momentum per unit mass 
flow. This quantity, which produces the useful torque, can be 


calculated by using the procedure described in section 3.5. 
Thus force due to momentum change = per unit mass 

flow. 


The mass flow per finite element 


m 

Ug. . Ug 


where m is the total mass flow and n^ is the number of finite 
elements (of equal length) in the blade, 


Therefore, force due to gas bending per element — 

(in Y-direction) 


(5.23) 
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The length of each finite elenent is given by — . figure 5.4 

H ^ 

shoT^s the degrees of freedom associated with a blade element 

and the nodal distances from the axis of rotation along with 

the end dimensions of the equivalent beam. The distance of the 
i/li 

i node from the axis of rotation (z^^) is given by 
1 

^i+1 = ^ 2^^ ~ ^ ^ ^ 5 i = 0,1,2,.. .,ng (5.24) 

H/ 


The whirl component of velocity varies with the radius and in 
such cases the free vortex relation is assumed to be valid so 
that 


(G 


V/ 


2 


-f c 


w. 



constant 


(5.25) 


By assuming that the end nodes share half the load 
as compared to the middle nodes, we can specify the gas bending 
loads corresponding to the i^^ degree of freedom by the 
following expressions s ■; ' 


P 


1 





+ 0 ) 

ruT * 


w. 


m 


4 )(. .,, 4 ^ ) 

^2 2^1 °E bj 


(5.26) 


j m /dw 

(n^n^+1) S “e “b ^ “e ^ 


-) (5.27) 


* / . ^ (n 4/^ >_S— ^ s i= ) ( 5 • 28 ) 
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where is the number of degrees of freedom at a node, v/hich 
is equal to 4 in the present case. Since all these forces 

s ^ in Y direction, they have to bo 

ili D 

resolved into components parallel to X and Y axes before 
solving the static equilibrium problem. 

(ii) lorces due to gas pressure: 

The total pressure force acting at the annulus will 
be (p 2 - p^ ) •Jrd h^. Thus the pressure force per element will 
be approximately 

(Po - p,) ird tu 

£ 2 ii fs qc 


Here again, we assume that the end nodes share half the load 
as compared to the middle nodes. We can specify the loading 
due to pressure at the end and middle nodes by the following 
relations; 


p - I ~ 

3-2^2^ 

_ p f (ng) tyl) ' hR^Pg - 

~ 2 ^ 2 ^ * ng n^ 


(5.50) 

(5.31) 


(i.n^+3) 


= TT Zj 


hjj(p. 


Ill 


“e “r 


i = 1 ,2. 




(5.52) 


Since all these loads P 3 , ^^,+ 3 ) ^ direction, 

"tkiey haYe "to be resolYed iuto componenljs parallel to X and Y 
3^es before solYing tb© static eg.uilibriuiE problem# 
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Once the load vector is computed, v/e can get the 
nodal deflc!Ctions corresponding to the various degrees of 
freedom of tho boam (in X and Y directions) by solving the 
equilibrium equations iK] U = P where [K3 is the assembled 


stiffness matrix, U is the vector of nodal displacements and 
P is vector of nodal loads. The maximum deflection will occur 
at the free end of the cantilever and its magnitude can be 


found as 


A 


t 


(u 


2 


+ u 


2 

(ngnp+3 ) 



(5.33) 


where u^^ are the displacement degrees of freedom of the 
structure . 

In Chapter 4, the displacement due to bending ) 
was expressed in terms of the vector of nodal displacements 
(Ui) of the element. For the first element of the blade in 
which the stress would be maximum, the bending deflection 
in the YZ piano is given by (by neglecting the effect of shear 
deformation) : 


w, (z) = -1 (2z^ - 31z^ + 1^) + “I (31z^ - 2z^) - 

^ 1 ^ 1 ^ 

^ (z^ - 21z^ + l^z) - -| (z^ - Iz^) (5.34) 

1 ^ 1 

where 1 is the length of the element, which is equal to 
hj^/ug and z is the coordinate measured along the Z-axis 
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(Figure 5.4). The stress and the strain e are related 

zz 


with as 


0 = E e 

zz zz 


yB 


2 

-iJ!h 


u. 


Sy {-4 (i2z - 61) 
1^ 


+ 


u 


3 (61 -- 12z) - (6z - 41) - (6z ~ 21)} 

1 1 “^ 1 ^ 


(5.35) 


whore y is the distance of the fibre tneasured along Y-axis 
from the neutral axis. 

The iuaxiiDum value of stress ( a ) will occur at 

zz max 

^ "" ^"^^eq^root/^ ^ = ^ is given by 


(cr ) 
zz 


zz max 
6u. 6Uc 4u^ 2u, 


max 


( 5 . 36 ) 


Since node 1 is fixed, u.j = U 2 = 0 and hence the equation 


( 5 . 36 ) simplifies to 


(CT ) 
zz 


max 




(hq) 


(5.37) 

By proceeding in a similar manner, the maximum stress induced 
due to the displacement of the beam in XZ plane can ba obtained 
as 

E (b_) 


(n ) 

ZZ'^ 


max 


1 


root + ug) 


{5.38) 
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Thus the total maximum stress induced due to gas bending and 
pressure force (a^) can be found as 

'O 


E(t ) 
eq 


root 


1 


3Uf- 

+ 


E(b 


u, 


;)! 


eq' 


root 


1 


3u„ 


+ Uq)| (5.39) 


Finally the total maximum stress induced at the root ( a ) 

inax^ 

will be 


max 




(5.40) 


max 


where is give n^ by equation (5.39) and (cr^^) hy equation 

( 5 . 20 ) . 


max 


5.3 Deflection, Stress, Vibration and Other Constraints 

By idealizing the blade as a rectangular, tapered, 
twisted and rotating cantilever beam, its deflection, stress 
and the fundamental natural frequency can be found using the 
finite element analysis as described in sections 5.1, 5.2, 4.4 
and 4 . 5 . Then the mechanical constraints can be evaluated in 
the manner indicated in section 2.5. The other behavior 
constraints of section 2.5 can also be evaluated by using 
the values computed in Chapter 3. 
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CHAPTER 6 

OPTIMIZATION METHOD 

In Chapter 2, the problem of optimum design of gas 
turbine stage has been formulated as a nonlinear mathematical 
programming problem. The methods of computing objective function 
and constraints have been developed in Chapters 3, 4- and 5. Ibe 
choice of the method of optimization and its description is the 
topic of this chapter. 

6.1 Choice of the Method 

The three general classes of widely used nonlinear 
programming methods are as follows; 

(a) Gradient projection method of Rosen^^ which was subsequently 
modified by Goldfarb^^ Though this method works well with 
linear constraints, its efficiency is considerably reduced 
in the case of nonlinear constraints. 

(b) Eeasible direction method of Zontendijk^ : This metljod is 
based on the generation of usable feasible directions at 
constraint boundaries. Although this method works in a 
direct manner in solving the problm, the analyses during 
optimization have to be done accurately as they influence 
the rate of convergence and accuracy. 
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(c) Penalty function methods; These methods are quite reliable 
and their sequential nature allows a gradual approach to 
criticalitj’- of constraints. The methods allow coarse- 
approximations to be used during early stages of tte opti- 
mization procedure and finer or more accurate approximations 
during the final stages. The computational time, however, 
is expected to be slightly more in these methods as the 
minimization problem has to be solved a number of times. 

If computational time is to be reduced, the penalty function 
methods allow the use of approximate analysis without 
effecting the accuracy and stability of the procedure. In 
the present work the penalty function method of Piacco and 
McCormick^"^ is used as it has. been found to be quite 
reliable • 


6.2 Penalty Function Methods 

Penalty function methods transform the basic problem 
into alternative formulations such that numerical solutions are 
sought by solving a sequence of unconstrained minimization 
problems . 

Let the basic optimization problem be of the form: 


Pind X such that 

f(X) minimum 

and 




( 6 . 1 ) 


3 = ■> 5 2 , 


• • . j 


in 
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Thio ppoblBin is con.vei’tisd into an unconstrainsd ininiiiiiza'tion 
problem by constructing a function of the form 

IQ, ^ 

0(X, r^) = f(X) + ^ a g (X) (6.2) 

0=1 

where G- is some function of the constraints. 

If the minimization of the ^-function is repeated for a sequence 
of values of the response factor r^^, the solution may be brought 
to converge to that'^of the original constrained problem given by 
equation (6.1). 


6.3 Piac co-McCormick Interior Penalty Function Method 

In this method, the objective function is augmented 
with a penalty term consisting of the constraints as 


0(1 r,) = f(X)-r, ^ (6.3) 

t] 

Where 0(X, r^) is called the penalty function. The minimization 
of jZf-function is performed for a decreasing sequence of r^ so 
that 

r < r. (6.4) 

^k+1 ^ k 

Equation (6.3) requires a feasible starting point and, in the 
present work, it is found by a process of trial and error. 

Since each of the designs generated during the optimization 
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process lies inside the acceptable design space, the method is 
classified as interior penalty function formulation. 


6.4 David on-Pletcher-Powell Variable Metric Unconstrained 
Minimization Method 


In the penalty function formulations, since a sequence 
of unconstrained minimizations has to be performed, the selection 
of an efficient method of unconstrained minimization becomes very 
important. All the methods of unconstrained minimization find a 
sequence of improved approximations to the optimum according to 
the iteration; 


‘i+1 


X. + T s. 
1 1 


(6.5) 


where 

X. . = the design vector corresponding to the minimum of 

^-function along the current search direction 

X. = the starting design vector 
^ -*■ 

T* = the minimizing step length in the direction S^. 

There are several methods available for finding the search 
direction in equation (6.5)- 

In the present work, the David on-Pletcher-Powell 
Variable Metric method^® is used. This method can be considered 
as a quasi-Uewton algorithm, and is a very powerful general 

local unconstrained minimum of a 


procedure for finding a 
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function of many variables^^. In this method, the i"^^ search 
vector in equation (6.5) is computed as follov^si 


S. 


[H, ] V^(X) 


( 6 . 6 ) 


where denotes the gradient of the jZl-function at and 
[H^J is a positive definite symmetric matrix. The matrix Ch^] 
is updated according to the following procedure: 


and 


= [Hj,] + Ca.] + [B.] 


(6.7) 


where 


[A,] 


T S. S. 
1 1 


(6 . 8 ) 


[B,] 


[ H^] f 


(6.9) 




-VJZ({X,) 


( 6 . 10 ) 


The updating of Ch^J preserves the symmetric positive 

definiteness of which ensures the stability of the 

procedure. The positive definiteness of is influenced 

only by the accuracy with which t is determined. To start 

with [H ] is taken as the identity matrix in the present work, 
o 
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5 5 One Dimensional Minimization Method 

In the present wo he, the one dimensional minimization 
is accomplished by cubic interpolation method. In this method, 
the algorithm used to compute t* is reapplied until the cosine 
of the angle between the vectors S. and at the minimising 

step length t* is sufficiently small (e), i.e. 


. V )Z( 


COS© = 


iil. 


( 6 . 11 ) 


S. 

i' 


V jZS, 


i+r 


Y is the minimum along the 
This ensures that is rut 

In the present work, the value of e is taken as 
to reduce the computer time, the number of cubic 


direction S^. 
0.05. However 
int e rpol at ions 


is limited to three. 


6.6 Additional Considerations and Convergence Criteria 
(i) Starting point 

r.-e 0^(T r ), the starting feasible 
For the minimization of j2>(X, 

point X„ is found by a process of trial and error. Each 
subsequent stage used the solution of the previous sta^ es^a 

starting point. In some oases, the overall procedure 

oY+-rar)olation technique to 
accelerated by employing an extrapola 

•.-H. for subsequent unconstrained minimi 

determine starting points have been 

^ m-i nimization stages have oeen 

satlon cycles (after two or mo „oXation mat be 

completed ) . Starting points ob 



141 


checked to ensure that they satisfy the constraints, because at 
each st.ige, it is necessary to start the unconstrained minimi- 
zation of 0(X, from an acceptable design point. 

(ii) Initial value of r^ 

If r^ is largo, the function is easy to minimize, but 
the minimum may lie far from the desired solution of the 
original constrained minimization problem. On the other hand, 
if r^ is small the function will be hard to minimize. In the 
present work, the value of r^ is chosen such that 

1.25f(X) 1 0(^0’ 1 2.00f(X^) (6.12) 

(iii) Subsequent values of r^ 

The total number of to be employed is given as an 

input to the problem and the values of found by using 

the ratio 

^^+-1 0.1 (6.13) 

^k 

(iv) Initial positive definite matrix ] 

The identity matrix [I] is used for the initial [H^]* 

(v) Updating the, [H^J matrix 

Whenever the cosine of the angle between S and is 
less than -0.1, the matrix is taken as [H^J. Otherwise 

equation (6.7) is used to find [ • 
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(vi) Restarting the [h] matrix 


In the case of highly distorted or eccentric functions, 
it might occur after few iterations that S? . v0. is positive, 
indicating that is not a direction of descent. When this 
happens, the remedy is to set [h^] back to and proceed as 

if starting over again. 


(vii) Termination of minimization of each r^ 

For each r^, the minimization of the JZ^-function is 
terminated whenever the predicted percentage difference between 
the current and the optimal 0-values is less than a small ntunber. 
Thus the convergence criterion can be expressed as 


[H.] V0. 


S^,0. 

1 


(6.14) 


where e is a small quantity. The value of e used in the present 
work va.ded from about 0.10 for r^ to 0.0005 for higher r^^. 


(viii) Gradient of 0-function 

The gradient of the 0-function can be expressed as 

V0 = vf + r I ^ . Vgj_ 

0=1 gj 

However the gradient of the 0-function has been calculated by 
using the backward difference formula in the present woik. An 
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efficient method of evaluating the gradients of static deflection, 
natural frequencies and eigen vectors has been suggested by 
Eeddy 


(ix) Number of cubic interpolations 


9 

A maximum of five cubic interpolations are allowed in 
each one dimensional search problem. Out of these, only the 
final interpolation involves the evaluation of the gradient ^0. 
All preliminary interpolations use a perturbation scheme to 
determine the dot product S ..vP as 

gT rf ^ - U f i (6.16) 

{ t"" - D 


where ani define the step length. 


(x) Kuhn-Tucker conditions 


Before stopping the minimization process of ary problem, 
the Kuhn-Tucker conditions are tested as the necessary conditions 
for the optimum. Mathematically these conditions can be expres- 


sed by; 



+ I 

je J 


0 , i = 1 , 2 , 


X. > 0 
3 " 


3 e J 


(6.17) 

(6.18) 


where X. are called the Kuhn-Tucker multipliers and J is the 
d 

set of indices of active constraints. It is important to not 
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that the Kuhn— Tuck©r conditions are necessary and sufficient 
for a global niininiuni only in the case of convex progrannntng 
problems. However, equations (6.1?) and (6.18) can be used as 
the necessary conditions to test the minimum of any practical 
design problem. 

(xi) Relative minima 

In order to see whether any relative minima exists in 
the design space, two completely different starting points may 
be used for the sequence of minimizations for ary example. If 
the two sequences lead to the same optimum design (except for a 
small difference that might occur due to numerical instability), 
it can be assumed that local optimum is same as the global 
optimum. 

(xii) Reducing the total computational time 

It has been observed that seme of the automated 
optimum design problems take longer tine to satisfy the prescri- 
bed convergence criteria even after reaching very near to the 
optimum design point. This happens whenever the function being 
minimized is highly distorted or eccentric. In such cases, it 
may not be worthwhile to try to reach the exact minimum to 
obtain about 0.5 or decrease in the objective at the expense 
of 40 to 50^ more computing time. This type of situation can 
be identified by manually inspecting the progress of the 
optimization path at various stages. 
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CHAPTER 7 

iroijiERIOAIj EXMPLES MD SEISITIVITI AHAPYSIS 

A computer program has been developed to implement 
the optimization procedure described in Chapter 6. The program 
contains twenty one subroutines along with the main program. 

The details of the subroutines and their calling sequence are 
discussed in Appendix E. The program developed is quite general 
and can be used for the optimization of ary axial flow turbine. 
Other problems in the area of design of turbomachines can also 
be solved by making slight modifications to the present program. 
The flow diagram for Cholesky decomposition of symmetric banded 
matrices used to obtain a partial solution to the eigen value 
problem has been also shown in Appendix E. To demonstrate the 
effectiveness and flexibility of the program developed, the 
design of an axial flow gas turbine stage is considered. The 
problem is solved with different objective functions under 
different constraint sets. The data pertinent to the design of 
the stage along with the material properties of blades and gas 
properties are given in Table 7.1. 

7.1 Maximization of Isentropic Stage Efficiency 

The maximization of isentropic stage efficiency or 
minimization of losses is considered as the first example. The 
effects of rotation, taper and twist of the blade are considered 
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in the finite element idealization. Since the effects of shear 
deformation and rotary inertia on the first natural frequency 
are expected to be very small, they hare been neglected during 
computations to secve computer time (by reducing the nodal 
degrees of freedom by fifty percent). The components of the 
design vector are taken as: 


= Mean diameter of rotor 

Xg = Ratio of the chord of rotor blade to mean diameter 
X, = Ratio of the chord of nozzle blade to mean diameter 
X^ = Ratio of spacing to diameter, at mean radius of the nozzle 
blades 

X^ = Ratio of spacing to diameter at mean radius of rotor blades 
Xg = Relative angle of the velocity triangle at inlet of the 
rotor blade at mean radius 

Xy = Relative angle of the velocity triangle at the exit of the 

rotor blade at mean radius 
Xg = Axial velocity of flow across the stage. 

The constraints of the problem can be stated as follows. 
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(7.2) 
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(7.28) 

(7.29) 

(7.30) 


The 

X, 


-p +Vir> desisn variables are tahoii S'® 
starting values of the ^ ^ 

= 0.432, = 0.046, = 0*0^°’ 4" * ’ ^ 

= 20.5°, = 54.57° and Xg - 272.0. 
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The upper and lower bounds as well as the optimum 
values of the design variables are shown in Table 7.2. The 
bounds on the behavior constraints and the values of the 
response quantities at the starting and optimum designs are 
shown in Table 7.3. The progress of the optimization path, 
showing the cumulative number of one dimensional minimizations 
versus objective function (losses or one minus efficiency) is 
shown in Figure 7.1. In this figure, the variation of penalty 
function and mass of the stage has also been represented. The 
results of optimization show that there is a 30.5% reduction 
in the objective function. The optimum point corresponds to 
an increase in the efficiency of the stage by 2.48% and a 
reduction in weight of the stage by 55.3% compared to the 
starting design vector. None of the side constraints is active 
at the optimum point. The constraints on stress and degree of 
reaction at root have become active (out of the behavior 
constraints) at the optimum point. It can be obseiwed that 
X- , X, and X^ have not changed much from starting point while 
the other variables show appreciable change from the initial 
starting point. The optimum point has been found with 3 values 
of r^ in 15 one dimensional minimization steps which required 
about 61 minutes of computer time on an IBM 7044 compiter. 

Although the point obtained at the end of 15 one 
dimensional steps (with 3 values, of r^^) can be taken as the 
optimum (evident from Figure 7.1 )> minimization is carried 
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for throe more values of r^^ which required an additional 
computational time of 70 minutes. The results obtained with 
6 values of are compared with those obtained with 5 values 
of r^^ in Tables 7.2 and 7.3. It can be seen that while the 
first 3 values of r^^ reduced the objective function (one minus 
efficiency) from 0.0813 to 0.0565, the additional 3 values of 
r^ could reduce it to 0.0557 only. Thus no significant reduction 
in the objective could be achieved by proceeding beyond 3 values 
of r^^. Since this behavior is characteristic of the interior 
penalty -function method, all the subsequent examples are 
solved by using 3 values of r^ only. 

In order to test whether the optimum point found 

corresponds to a local minimum or the absolute minimum in the 

design space, the same example has been solved with a second 

starting design vector whose components are = 0.276, = 

0.056, = 0.068, = 0.059, X^ = 0.028, Xg = 18.6°, X^ = 

52.5° and X^ = 252.7. The progress of the optimization path 
o 

is shown in figure 7.2. This plot is similar to figure 7.1 and 
the optimum point obtained (f = 0.0543) also compares well 
with that of the previous case (f = 0.0557). The starting and 
the optimum designs of the two cases are compared in Tables 
7.2 and 7.3. The optimimi design variables are also xn good 
agreement with each other. Apart from a small difference that 
might have occurred due to numerical instability, the two 
optimum points appear to be same. Although, merely on the 



151 


basis of two trial starting designs, it is hard to say that the 
minimum obtained is the absolute minimum over the design space; 
finding the similar optimum design by starting from t7;o diff- 
erent initial designs is at least a pointer in that direction. 


7.2 Minimization of Weight 


In the second example , the optimization of weight 
(mass) of the axial flow gas turbine stage is considered. The 
effects of rotation, taper and twist of the blade are considered 
In the deflection and vibration analysis of blades. The cons- 
traints, design variables and their bounds are same as in the 

case of the first example (first starting point). 

The optimisation results are shown in Tables 7.2 and 

7.3. The progress of the optimisation path, show'ng 
amative number of one dimensional minimisations versus the 

Objective function is shown in .l^re 7.3. The optimisation^ 

n -o o ffn reduction in the objective, 

results show that there is a • ^ 

, i nut of all the side and behavior constraints, 

function (mass), ait of a 

only the degree of reaction at roo 

ri-t-Tirp t;h8 poxnti* 

of the rotor blade have become a their 

„ V ud l' have not changed much from their 
Comparatively Xg , and g 7 , it can be observed 

initial starting values. Jrom Pigure 7.3 it 

e Uv email a»unt first and then 
that the losses decrease y 

increased at a faster ra e. efficiency. Por 

in weight can be obtained only a 
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■this sXciinplG j "th© conip'utG]? tini© rsQ^uirsd, is 60 niinutss on an 
IBM 7044 computer for 15 one dimensional minimization steps. 

7.3 Optimization of Weighted Combination of Efficiency and Mass 

The minimization of a weighted sum of efficiency and 
weight (mass) of the gas turbine stage is considered in the 
third example . The constraints, design variables and their 
bounds are same as in the case of example 2. The blade has 
been idealized as a rotating, tapered and twisted cantilever 
beam by neglecting the effects of shear deformation and 
rotary inertia. The magnitude of losses and mass (weight) of 
the stage are normalized such that their contributions will be 
equal at the starting design point. The optimization results 
are given in Tables 7.4 and 7.5. The progress of the opti- 
mization path is shown in Eigure 7.4, where the variations of 
penalty function, losses and mass of the stage with the number 
of one dimensional minimizations are shown. It is observed 
that there is a 45?^ reduction in the objective function while 
the efficiency increased by 1 .65?^ and the weight reduced by 
69.5%. The variables , Xg and have not changed much 
from the starting point. The constraints on the degree of 
reaction at the root and the stress have become active at the 
optimum point. It can be seen that the reduction in losses 
or mass of the stage is less in this case compared to the 
case where the individual objective function xs considered 
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for optimization. Thus the present results indicate a compro- 
mise between efficiency and mass of the stage. Here 15 one 
dimensional minimization steps recuired about 62 minutes of 
computer time on an IBM 7044 computer. 

7,4 Optimization of Weighted Combination of Mass and Efficiency 
(Without Considering Rotation) 

In the fourth example , the minimization of vfeighted 
sum of efficiency and mass (weight) of the stage is considered. 
The rotor blade has been idealized by using four tapered and 
twisted beam elements. The effect of rotation has not been 
included while computing eigen values of the idealized blade. 

The shear deformation and rotary inertia effects have also been 
neglected . 

The optimization results of this problem are shown 
in Tables 7.4 and 7.5. The progress of the optimization path, 
shovi^ing the values of penalty function, objective function, 
losses and weight versus the number of one dimensional minimi- 
zation steps is shown in figure 7.5. In this case, the 
objective function has decreased by 43^ > efficiency has 
increased by 1 .68% and the mass of the stage has decreased by 
S5.5%. Most of the variables behaved as in the case of third 
example. The degree of reaction at the root and the stress 
at the root have become active at the optimum point. The 
program required about 61 minutes of computer time for 15 one 
dimensional minimization steps. 
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7.5 Optimization of Weighted Combination of Efficiency and 
Mass by Considering Only the Side Constraints 

In the fifth problem, a linear combination of the 
mass and the losses of the stage has been minimized with equal 
contribution of losses and mass at the starting design point. 
Only the side constraints (g^ to g-|g) are considered with no 
constraints on the bohaTior (response) quantities in the 
problem formulation and solution. The blade has been idealized 
with four finite elements and the effects of rotation, shear 
deformation and rotary inertia have been neglected. 

The optimization results are shown in Tables ,7.4 and 
7.5. The progress of optimization path, showing the cumulative 
number of one dimensional minimizations versus the values of 
objective function, losses and mass is shown in Figure 7.6* 

It is found that the objective function was reduced by 47.6^ 
while there is an increase of 2.1?^ in efficiency and a 
reduction of 69.5?^ in mass of the stage. It is observed that 
although the objective function (losses as well as mass of the 
stage) has been reduced by a larger amount, the behavior 
quantities like stress, degree of reaction at root and tip 
deflection have become larger than the permissible values. 

This clearly demonstrates the necessity of putting bounds on 
these response quantities. The program required only ^out 
20 minutes of computer time on IBM 7044 computer for 15 one 
dimensional minimization steps in this case. 
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7.6 Sensitivity Analysis 

In practice a designer would be interested in knowing 
how the response quantities vary with a change in the design 
variables . This type of sensitivity analysis will help the 
designer in manipulating the design variables to suit some 
specific requirements. Further, in some cases, the results 
obtained from the optimization procedure may have to be rounded- 
off to the nearest practical values of the design variables. 
Hence a sensitivity analysis of response quantities, namely, 
objective function, losses, mass, stress, degree of reaction 
at root and moan radius, first natural frequency, flow coeffi- 
cient, stage temperature drop coefficient and angle of nozzle 
blade at outlet, with respect to the various design variables 
is conducted. In this analysis, the reference design is 
taken same as the optimum point of example three (section 7.3)* 
The design variables are varied on the negative and positive 
sides of the reference (optimum) values and the magnitudes of 
the behavior quantities are plotted against the percentage 
charges of the design variables in Figures 7.7 to 7.15. 

From Figure 7.7 it is observed that the objective 
function (weighted combination of losses aod mass) is quite 
sensitive to variations inX^, Xgj andX^. The sensitivity 
of the two components of objective function is shown in 
Figures 7.8 and 7.9. Figure 7.8 shows that losses are more 
sensitive to the variables Xg > X^, X^ and X^ compared to the 
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other four variables. Similarly Figure 7.9 indicates that mass 
of the stage is more sensitive to the diameter of the turbine 
compared to the variables and Xg. The mass is almost 

independent of the variables X^, X^, X^ and Xg. 

The sensitivity of mechanical response quantities, 
namely, stress, deflection and fundamental natural frequency 
with respect to the various design parameters is shown in 
Figures 7*10 to 7.12. It can be observed that all these qua.n- 
tities are most sensitive to the variables X^ and X^ and less 
sensitive to the other variables. 

Figures 7.13 to 7.15 represent the variation of the 
aerodynamic response quantities, namely, degree of reaction at 
root, flow coefficient, stage temperature drop coefficient, 
stator blade angle at outlet and degree of reaction at mean 
radius with a percentage change in design variables. It can 
be seen (Figure 7.13) that the degree of reaction at root is 
more sensitive to the variables X,^ and Xg compared to the 
other variables. The stage temperature drop coefficient is 
found (Figure 7.14) to be qufte sensitive to X^ , Xg, X^ and Xg, 
and independent of the variables Xgj X^ and X^. The flow 
coefficient is found to increase (decrease) with increasing 
values of Xg(X^) while it is independent of the remaining six 
variables. Figure 7.15 represents the sensitivity analysis of 
degree of reaction at mean radius and gas outlet angle of 
stator blades. It is seen that the degree of reaction is most 
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sensitive to the variable X , less sensitive to X^, X. and X^, 

I loo 

and insensitive to the variables X 2 , X^, X^ and X^. Similarly 
the stator blade outlet angle (a^) is found to be more sensi- 
tive to X^ , Xg and Xg and insensitive to the remaining -variables 
(Xg > j ’ ^5 and X^) , 

Similar results have been obtained when the optimum 
point of example 1 (first starting point) is pertin?bed. These 
results are shown in tabular form in Table 7.6. 
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Tmm 7.1 

Data for Optimization Problems 

Number of design -variables, n = 8 
Number of constraints = 30 
Number of finite elements in a blade = 4 
Number of eigen values computed = 2 

Number of degree's of freedom of an element = 8 (since shear 
deformation effect is neglected) 

Maximum number of cubic interpolations used in one dimensional 
minimization = 3 

Maximum number of unconstrained minimuzation steps for any value 
of penalty parameter = 5 

Gas constant, R = 287 N-m/Kg-K 

Viscosity, "y = 4.0 x 10 "^ Kg/m-sec 

Ratio of specific heats of the gas, y = 1.333 

Density of the blade material, = 800 Kg/m 

11 2 

Young's modulus of the material, E == 2.0 x 10 N/m 
Shear modulus, G = B/ 2,6 

Stagnation pressure of enterir^ gas in the stage, 

= 4.0 X 10^ N/m^ 

Stagnation temperat-ure of entering gas in the stage, 

= 1100 K 

Mass flow of the gases across the stage = 20 Kg/sec 
Rotational speed of the turbine, N = 250 rps 
Radial clearance of the blades, k = .000154 ni 
Breadth taper ratio, 3 = 2.0 
Depth taper ratio, o = 1 *5 
Shear coefficient , v = .833 
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TABLE 7.3 

Rdsponse Quantities at Initial and Optinun Points for Example 1 and 2 
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Initial and Optinun Eesponse Quantities for Minimization of Weighted 
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CHAPTER 8 

CONCLUSIOHS AID EEC OMffiUDAT IONS 

lh .0 i'03,taibili'ty of au 1)011131; 6 d optiniuin dosign of axial 
flow gas turbine stage has been demonstrated by the numerical 
examples oi Chap cor ] , In this chapter conclusions drawn from 
the present study arc stated and the scope of extension of the 
present work is suggested. 

8.1 Conclusions 

(1) Computation of the efficiency of gas turbine stage 

The convergence of the iterative method developed for 
the calculation of efficiency has been found to be very good. 

The method converged within 3 to 4 iterations irrespective of 
the initial trial value given. With higher convergence 
criteria, the number of iterations may increase by one or two. 

(2) Computorization of air property relations 

The polynomial equations obtained for evaluating the 
properties of air aro quite accurate and the maximiM error 
involved in most of the equations does not exceed 0.0535 percent. 
These equations would bo useful in ali analyses involving gas 


flow. 
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(3) Vibration and stress analysis of rotor blades 

The turbine blade of ai2?foil section has been conv- j. 
into an eq.uivalent doubly tapered and twisted rectangular 
cantilever beam for the purpose of analyzing deflection <,+ 

’ ^cregg 

and natural frequencies of vibration of blade. A new finito 
element has been developed for the stress and vibration 
analysis of tapered, twisted and rotating beams by including 
the effects of shear deformation and rotary inertia. The 
element developed has been found to give reasonably good results 
even w’ith four elements. 

The results of the tapered beam with varying depth and 

breadth taper ratios were quite comparable Yirith those reported 

28 

by Mabie and Rogers . The results of twisted uniform beams 

35 

were in good agreement v^fith those given by Rosard . The 
following observations have been made regarding the vibration 
characteristics of cantilever beams: 

(a) The effects of shear deformation and rotary inertia reduce 
the natural frequencies of vibration of beams. 

(b) The natural frequencies increase with an increase in the 
offset (Figure 4.6). 

(c) In certain mode shapes, the natural frequencies increase 
while in some cases, they decrease with an increase 

in the rotational speed of the beam. 

(d) An increase or decrease of depth or breadth taper ratio 
changes the frequency ratio of the beam. 
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(e) The tvN^ist also effucts in two different v/ays; it may 
increase or decrease the frequency ratio of the beam 
depending on the amount of twist and the mode shape. 

(f) The effects of shear deformation, twist, offset and taper 
ratio have been found to be more predominant on higher 
modes of vibration. 

The Cholesky decomposition of symmetric banded matrices 
has been used efficiently in solving equilibrium equations and 
in obtaining a partial solution to eigen value problem by 
Raleigh-Eifz subspace iteration algorithm. 

The stresses due to pressure force, gas bending and 
centrifugal action have been calculated individually. At the 
root of the blade , all the three stresses are added to get the 
maxinmi stress. The stresses due to pressure and gas bending 
were calculated by using the finite element analysis. The 
centrifugal stresses were calculated by using the conventional 
method. The orders of magnitude of these stresses have been 
found to be same as the ones reported in reference 93 for a 
specific case. 

( 4 ) Results of optimization 

The problem of design of axial flow gas turbine stage 
has been formulated as a nonlinear programming problem by 
including the constraints due to mechanical and aerodynamic 
considerations. The interior penalty function approach with 
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tho root and the stresses in the blade have been found to be 
active at the optimum point . 

By considering the savings realized through the optimi- 
zation procedure (between 1.65^ to 2.48?^ increase in efficiency 
of the stage and 65.5?^ to 80.5/^ reduction in weight in the 
examples considered), the computer time of about 60 minutes 
required for solving one optimization problem is justified. 

Thus the present procedure gives the designer complete informatioi 
automatically in the initial design phase. 

(5) Results of sensitivity analysis 

The results of sensitivity analysis are expected to 
be useful in rounding off the optimum design variables to the 
nearest available values. These are also useful in determining 
the relative importance of the various design variables. This 
information would be helpful in eliminating some of the less 
important design variables in future design studies. It has 
been observed that the losses are most sensitive to the chord 
and spacing of rotor blades while the mass of the stage is 
most sensitive to the mean diameter of the rotor. The 
mechanical response parameters, namely the deflection, stress 
and fundamental natural frequency of vibration of the blade, 
have been found to be most sensitive to variations in mean 
diameter of the rotor and the exit an^e of the rotor blade 
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(Bj). The degree of reaction at root was found to be most 
sensitive with respect to and axial velocitj of flow. 

(6) Computer program 

The computer program developed for the automated 
optimum design of axial flow gas turbine stages involved 21 
subroutines. It is quite general and can be used for solving 
any axial flaw gas turbine design problem with different 
objectives and constraint sets. 

8.2 Recommendations for future V/ork 

1. The loss due to trailing edge thickness, entrance losses, 
disc friction and windage losses, mechanical losses a nd 
partial admission losses (in the case of turbine with 
partial admission) can also be considered and the overall 
efficiency of the turbine can be considered for optimi- 
zation. However, since the present stage of knowledge does 
not permit quantification of these losses in terms of the 
design variables of the problem, more research need to be 
done in this direction. 

2 . To include the three dimensional effects in the analysis of 
turbine stage, other methods like stream line curvature 
method, actuator disc theory etc. can be used instead of 
the simple free vortex theory used in the present work. 
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3, The axial displacement degrees of freedom can also be 
included in the finite element formulation to find the 
centrifugal stresses induced in the rotor blade more 
accurately. The other types of stresses like thermal 
stresses, stress due to root fixing and stresses due to 
creep and fatigue can also be included in the analysis. 

4, The torsional nodal degrees of freedom can also be 
considered in the finite element analysis to find the 
coupled bending -bending -torsional frequencies of vibration 
of turbine blades. 

5* No experimental results are available on the natural 

frequencies of vibration of beams including the effects 
of rotation, taper and pretwist .■ Some experiments could 
bo conducted to predict the vibration response of such 
beams . 

6. The vibration analysis of the turbine stage can be made 
by considering the assembly of rotor disc and blades 
with root fixing and shrouding. 

7. The present optimization method can be extended to find 
the optimum design of multistage gas turbines. 

8. The relative efficiencies of optimization can be studied 
by solving the problem of optimum design of gas turbine 
stage by using other optimization methods like gradient 
projection method, Zoutendijk's feasible direction method 
and sequential linear programming method. 
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9. The design problems can be formulated and solved by 
using a reliability framework. 

10. Some experimental test rig for gas turbines can be 

set-up and the validity of the optimum results can be 


studied. 
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APPEIJDIX A 

POLYNOMIAL EXPRESSIONS EOE PROPILE LOSS COBPPICIENTS 

USED IN SECTION 3-5*4 


The profile loss coefficient Yp(i,j) is expressed as 
a polynoinial function of the spacing to chord ratio, x, as 
follows; 

4 

Y^(i,j) = I c x^ (A.1) 

P k=o ^ 

where i = 1 for nozzle blades with = 0 

= 2 for impulse blades with = 3^, 
and 0 = blade angle in degrees . 

The following equations were used for evaluating Y^ in the 
present work. 


Yp(1,80) 


- 1.0236312 x^ + 2.5686004 - 2.1936137 + 

0.74141858 X - 0.021356594 (A. 2) 


Yp(1,75) 


Y (1 ,70) 

Jr 


1.1784926 x^ - 3-2665715 x^ + 3-4585936 x^ - 
1.6462663 x + 0.34334468 (A. 3) 

- 0.26502239 x^ + 0.70910158 x^ - 0.49822962 x^ + 
0.0090238 X + 0.09513364 (A. 4) 


Yp(1,65) 


= 7.1051656 x^ 


24-956305 x^ + 32.491226 x^ 


18.53773 X + 3.9333972 


(A.5) 
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Yp(l,60) 

Yp(l,50) 

Yp(l, 40 ) 

Yp(2,70) 

Yp(2,65) 


6.9573477 - 24.391842 + 31.685563 x^ - 

18.048353 X + 3.8250498 (A. 6) 

7.1509773 x"^ - 24.834378 x^ + 31 .946036 x^ - 
18.029084 X + 3.7862075 (A. 7) 

6.9582213 x^ - 24.177738 x^ + 31.093230 x^ - 
17.537165 X + 3.6792156 (A. 8) 

0.46636044 x"^ - 1.6216202 x^ + 2.1681230 x^ - 
1 .1824825 X + 0.35962101 (A. 9) 

- 0.7677269 x^ + 1.7245656 x^ - 1.0102869 x^ - 
0.00304808 X + 0.20849764 (A. 10) 


Yp( 2 , 60 ) 


Yp(2,55) 


Yp(2,50) 


Yp(2,40) 


= 0.94616804 x^ - 2,8471011 x^ + 3.3114753 y? - 

1.726889 X + 0.44034477 (A. 11) 

= - 0.06295425 x^ - 0.02127604 y? + 0.50723384 x^ - 

0.59957593 x + 0.27957267 (A. 12) 

= - 0.32384541 x^ + 0.87772575 x^ - 0.56541410 x^ - 

0.09070201 X + 0.19223416 (A. 13) 

= - 0.04099546 x^ - 0.14282837 x^ + 0.72273100 x^ - 


0.78531104 X + 0.32239876 


(A.14) 
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Each of the equations (A. 2) to (A.I 4 ) was obtained from 
five data points taken from the profile loss curves. The data 
points used in obtaining the polynomials are given in Tables A.1 
and A. 2. 



Profile loss Data for Nozzle Blades ( 6, 


201 


o 

II 



CM o O UD 00 

MDtOCMi-T- 

o o o o o 

• • • • * 

o o o o o 

00 O T- 

o • • • • 

O O O r~ ^ 

^ CM OJ O 
M:) tOv CM CM CM 

o o o o o 
o o oo o 

rAVD CO O V“ 

• * « • • 

O O o V- T- 


VD VD VD 00 CM 
P> VO CM CM 

iM O O O O O 

o • • • • • 

o o o o o o 

VD 



KVVO 00 O ^ 

• • « * • 

O O O V- ^ 

CO CO Ovo CM 

o o o o o 

• « • • « 

o o o o o 

tnvo CO O T- 

* • • , • 

0 0 0'<-T- 

O O 00 O CO 
ir\Lr\Lr> 

o o o o o 

• • • • * 

o o o o o 


tACD CO O ■r- 

• • ■ • • 

O O O v~ T- 


00 O 00 o o 

mm ltn VD 
O O CO o o 

o o o o o 


O O O Ln 
"M- mvo 00 CT^ 

10*90 

o o o o o 


CM 00 VD 00 O 
vjommmo- 

o o o o o 

• • * • « 

o o o o o 


T- o o o o 

'M- IX\VD (T) 




O o CO OVD 
CTiVD 

-r- O O O O 


o o o o o 


CM o o o o 
tTiLDr-ovo 

• • • • • 

O O O O T~ 


^ mm o o 

CM CJ>!>-a) Ch 

o oo 
CD O CD o* <0 


mo O O O 
K^m^-cJ^o 

• * * • • 

o o o o ^ 


o o VO 
mo 00 c3^o 

O Or* 

• • • * • 

O O O O O 


mo O O O 
mmi>-crkO 

o o o Or- 


O 00 00 

O O CM 


o o o o o 


mo mo o 

mmvD 00 o 

• • • • • 

o o o o 


t-CM VD CD CM 
IXV CM T- CM Lf\ 


o o o o o 


O O o 

t<^Lr\vD CO O 

• • • • * 

O O O O 'T- 


O CO O CM CM 
VD mmvD cn 


o o o o o 


o moo o 
mmtr^co o 

• • « • • 

O O O O r- 


CM m*^m 



CD 


CM m-M" m 

















202 


APPEEDIX B 

AIR PE0PE3RTY EELATIOES 


If X and j denote, respectively, the independent and 

"th 

the dependent variables, the n degree polynomial relation 
between the two can be expressed as: 

y = y(x) = p^ + p^x + p^x^ + ... + . 

By transforming the independent variable as: x = 

c 

x^(x), the above polynomial relation cai be rewritten as: 

y = yU^) = + ... + a^x“ 

In the present work, a^^, b^, c^, d^, (i = 0, 1, 2 ... 
n) indicate the various polynomial coefficients in different 
temperature raages. The various air property relations are 
summarized below; 


(1) Temperature in terms of enthalpy (x = 0.01 x) 

c 

5 X 

y = I ^i^c ’ 100°R < T ^ 1000 R 

i=o 

where a^ = 0.6665893 E 00, a^ = 0.4189995 E 03, 

a2 = -0.3435174 E 01 , a^ = 0.4291952 E 01 , 

a^ = -0.2118072 E 01 , a^ = 0.2575576 E 00, 

4 

y = ^ b.x^ , 1000°R <_T < 3000°R 

i=o ^ ° 


0.4189995 E 03 = 0.4189995 x 10^ 


(B.la) 


(B.lb) 
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where = -0.4048060 E 02, = 0.4655783 E 03, 

b2 = -0.1588451 E 02, = 0.7941960 E 00, 

b^ = -0.1018477 E -01, 

5 . 

^ = I ’ 3000°R 1 T 1 6500°E 

i=o 

where = 0.1628055 E 03 , °i = 0.3799953 E 03, 
Cg = -0 .3107125 E 01 , = 0.5445013 E -01 

(2) Temperature in terms of relative pressure 

= In (10 x) 

6 

y = I aj^x^ , 100°E 1 T £ 1000°R 

1=0 

where a^ = 0.2543419 E 03, a^ = 0.7281138 E 02, 

ag = 0,1041926 E 02, = 0.1013960 E 01, 

= 0,7763628 E-01 , a^ = 0.2131233 E-02, 
ag = -0.4833005 E-03 , 

6 

7=1 > lOOO^R < T < 3000°R 

i=o — — 

where b^^ = 0.3367258 E 02, b^ = 0.1873639 E 03, 

bg = -0.1345989 E 01 , b^ = -0.2795454 E 01, 

b^ = 0.1234340 E 01 , b^ = -0.1190855 E 00, 


(B.lc) 


(B.2a) 


(B . 2b ) 


bg = 0.4245929 E-02, 



204 


y = y c.xj , 3000°E < T < 6500°R 

ilo ^ ° - - 

where = 0.5368253 E 03 , = -0.1870057 E 02, 

Cg = 0.2573325 E 02, = -0.1423321 E 00, 

= 0.5124743 E-01, = 0.3583317 E-02, 

(3) Temperature in terms of internal energy (x =0.01 x 
5 

y = I a.xj , 100°R < T ^1000°R 

i=o 

where a^ = 0.1058974 E 01 , a^ = 0.5870682 E 03, 
a2 = -0.7926372 E 01 , a^ = 0.1550686 E 02, 

a^ = -0.1135502 E 02, a^ = 0.2040550 E 01, 

5 

y = y h.xj , 1000°R < T < 3000°R 

x=o 

where h = -0.5067941 E 02, b, = 0.6665731 E 03, 

o ’1 ’ . 

b^ = -0.3422699 E 02, b^ = -0.2883410 E 00, 

b^ = 0.5878343 E 00, b^ = -0.4505020 E-01 , 

5 

y = I o.x^ , 3000°R < T < 6500°R 

i=o 

where c = 0.3549173 E 02, = 0.6147055 E 03? 

= -0.2961354 E 02, = 0.2515192 E 01, 

c. = -0.1171622 E 00, = 0.2249486 E-02 


(B.2c) 

) 

(B.3a) 

(B.3b) 

(B.3c) 
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(4) Temperature in terms of relative volume 


X. = In (100/x) 



4 





7 = 

.1 ’ 

1=0 

lOO^E < T < 1000 R 


(B.4a) 

where 

a = 0.6253480 E 

0 

03, 

a^ = 0.2485310 E 03, 



ag = 0.4707777 E 

02, 

= 0.5013478 E 01 , 



a^ = 0.2423038 B 

00, 










7 = 

y b.x , 

.L 1 0 ’ 

1=0 

1000 

°R < T <_ 3000°R 


(B.4b) 

where 

= 0.6352536 E 

03, 

b^ = 0.2254187 B 03, 



b2 = 0.6870292 E 

02, 

b^ = -0.3726730 E 

01 , 



b^ = 0.8656939 E 

y c . X J , 

i£o ^ ° 

00, 

°R ^ T ^ 6500°R 



7 = 

3000' 


(■B.4c) 

where 

Cq = 0.7072512 E 

03, 

C.J = 0.1925400 E 03, 



©2 = 0.6989422 E 

02, 

= -0.2417352 E 

01 , 



= 0.7004906 E 

00, 




(5) 

Temperature in terms of entropy (x^^ = e ) 



7 = 

? 1 
. l ’ 

1=0 

100° 

R < T < 1000°R 


(B.5a) 

where 

a = -0.5001787 E 01 , 

0 

= 0.8689595 E 

01 , 



ag = 0.3395702 E 

02, 

a, = 0.7477182 E 

01 , 
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= -0.1858278 E 01 , = 0.3150205 E 00, 

ag = -0.2599161 E-01 , 

6 

= I ’ 1000°Rj< T £ 3000°R (B.5b) 

1=0 ^ 

where = -0.1702268 E 03 , b^ = 0.8584996 E 02, 

bg = 0.2114853 E 02, b^ = 0.8462060 E 01, 

b^ = -0.1273866 S 01, b^ = 0.6457035 E-01, 

bg = -0.6959124 E-03, 

4 . . 

j = I , 3000°R <T < 6500°R (B.5c) 

1=0 

where = -0.1293062 E 04 , = 0.4627224 E 03, 

Cg = 0.1436791 E 01,, = 0.9997245 E 00, 

= -0.2324186 E-01 

(6) Enthalpy In terms of temperature (x = .01 x) 

4 . 

y = Z ^l^c ’ 100°R < T < 1000°E (B.6a) 


x=o 


where a^ = -0.1896111 E 00, 

a^ = 0.2392554 E 02, 

ag = 0.8177288 E-02, 

a^ = -0.3637402 E-02, 

= 0.4739440 E-03, 



3 

= y b.xj , 1000°E < T < 3000°E 

x=o 


y 


(B.6b) 
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where h^ = 0.1380925 E 02, h-, = 0.2042364 E 02, 

bg = 0.2533055 E 00, = -0.2364187 E-02 , 


y = .1 


1=0 


•^c ’ 


3000°R <_ 6500°R 


where = -0.3553479 E 02 , = 0.2542318 E 02, 

= 0.8214479 E-01 , = -0.3846921 E-03 

(7) Relative pressure in terms of temperature (x^ = ( 


V -y^ 

y = .1 ^i^c 


100°R <11 400®R 


L “i c ’ 
i=o 


= 0.1537417 E-02, 

a^ = -0.2998780 E-02, 

ap = 0.2099430 E-02, 

a, = -0.9920470 E-04, 
3 

= 0.5367652 E-03, 

a^ = -0.4252450 E-04, 
5 

a. = 0.1781536 E-05, 

6 



y = 


^ i 
1=0 


400°R <.1 ^1400 


where b^ = 0.1704886 E 00, 
bo = 0.6571773 E-02, 

C, 

b, = 0.9983916 E-04, 
4 

bg = 0.2577276 E-07, 


b^ = -0.6836948 E-01 , 
b^ = 0.1621208 E-02, 
bg = 0.1127265 E-05, 


5 

.1 

1=0 


y = .1. "i^c ’ 


1500°R ill 3000°R 


(B.6c) 

.015 x) 
(B.7a) 


(B.7b) 


(B.7c) 
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where = 0.2833894 E 01 , = 0.1159899 E 00, 

Cg = -0.3264126 E-01, = 0.3273487 E-02, 

= 0.6062104 E-04, = 0.-2453331 E-05 , 

6 

j = I d.xj , 3000°R < T < 6500°R (B.7d) 

i=o 

where = -0.1474357 E 03 , d^ = 0.1026554 E 02, 

dg = -0.2287566 B 00, d^ = 0.2089396 E-02, 

d^ = 0.1446234 E-03, d^ = 0.1600547 E-05 , 

dg = 0.1381830 E-08 

(8) Internal energy in terms of temperature (x = 0.05 x) 

c 

5 

^ = .1 1 1000°R (B.8a) 

i=o 

where a^ = -0.1879993 E 00, a^ = 0.3408698 E 01, 

ag = 0.1413006 E-02, a^ = -0.9960654 E-04, 

a^ = 0.2502025 E-05, = -0.1451709 E-07, 

6 

j = I b^x^ , 1000°R < I < 3000°R (B.8b) 

i=o 

where = 0.1340142 E 02 , b^ = 0.2697001 E 01, 

b2 = 0.1166697 E-01, b^ = -0.5392805 E-04, 

b^ = 0.3659680 E-06 , b^ = -0.1849987 E-08, 

bg = 0.3669679 E-11 

5 

= I c.xj , 3000°R < T < 6500°R 

i=o ^ 


J 


(B.8c) 
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where = 0.2527138 E 02, = 0.2539249 E 01, 

02 = 0.1174141 E-01, = -0.3073682 E-04, 

= 0.3878235 E-07, = -0.1526544 E-10 

(9) RelatiTe volume in terms of temperature (x = 80 

c 

4 

j = I a.xj , 100°R < T <_400°R 

i=o 

where a^ = 0.5282485 E 02, a^ = -0,8309882 E 00, 

a2 = 0.8322977 E-02, a^ = 0.1109702 E-04, 

a^ = -0.1839174 E-08, 

6 . 

y = I ’ 400°R 111 1300°R 

i=o 

where = 0.1825'952 E 01 , = -0.2004256 E 00, 

bg = 0.6083203 E-02, b^ = 0.8056071 E-05 , 

b^ = 0.3655806 E-07, b^ = 0.7191527 E-10, 

6 . 

y = I c.xj , 1300®R 111 3000°R 

i=o 

where = 0.4826241 E 00, = -0.4728824 E-01, 

C2 = 0.9295711 E-03, = 0.7207075 E-04, 

= -0.5315893 E-07, = -0.3531268 E-08, 

cg = 0.1579739 E-10 

s 

= T d.x^ , 3000°R < T < 6500°R 

x=o 


OOO/x) 

(B.9a) 

(B.9b) 

(B.9c) 


y 


(B.9d) 
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where = 0.1891114 E-01 , d^ = -0.3385984 E-02 , 

dg = -0.3200193 E-04, d^ = 0.5159768 E-04 , 

d^ = 0.5339627 E-06, d^ = 0.1317653 E-07, 

dg = -0.4221174 E-09 

(10) Entropy in terms of temperature 
= In ( 0 . 1 X ) 

4 ^ 

y = iL ’ 100°E < T < 500°E (B.lOa) 

Where = -0.29292898 E 00, a^ = 0.15769787 E 00, 
ag = 0.40449134 E-01, = -0.87987516 E-02, 

a^ = 0.70923293 E-03 , 

4 . 

7=1 ’ 500°E < T < 2000°E (B.lOb) 

i=o 

where = -0.67293304 E 00, = 0.41571832 E 00, 

b2 = -0.15920750 E-01, b^ = -0.528 2 9841 E-02, 
b^ = 0.80601767 E-03 

4 i 

7=1 , 2000 E < T < 6500°R (B.IOc) 

i=o 

where = 0.24348077 E 00, = -0.29914621 E-01, 

C 2 = 0.26780947 E-01, = 0.17085150 E-02, 


= -0.20158477 E-03 
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( 11 ) 

Y = 

where 

Y = 

where 

Y = 

where 

( 12 ) 

Y = 


Specific heat at constant pressure in terms of temperature 
(x^ = 0.1 x) 

6 ^ 

.1 ’ 100°R < T < 1000°E (B.lla) 

1=0 

% = 0.23956114 E 00, a^ = -0.77241925 E-04, 

Bg = 0.58395085 E-05 , a^ = -0.20386447 E-06 , 

= 0.34366968 E-08 , = -0.24666058 E-10, 

ag = 0.65225817 E-13, 

6 ^ 

I j 1000°R <_ T <.3000°R (B.llb) 

bo = 0.25720694 E 00, b^ = -0.89639173 E-03 , 

bg = 0.14075657 E-04, b^ = -0.80044782 E-07, 

b^ = 0.23861391 E-09, b^ = -0.38169716 E-12, 

bg = 0.26582135 E-15, 

5 

Jlq ^±^0 ’ 3000°R < T < 6000°R (B.llc) 

Co = 0.23471488 E 00-, c^ = 0.26026914 E-03, 

C 2 = -0.41468497 E-07, c^ = -0.10828877 E-08, 

= 0.19215360 E- 11 , c^ = -0.10247290 E-I 4 

Specific heat at constant volume in terms of temperature 
(x^ = 0.001 x) 

3 

I ’ 100°R < T < 500°R (B.12a) 

i=o 
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where = 0.17014188 E 00, = 0.91932141 E-02 , 

a2 = -0.42341176 E-01 , a^ = 0.57386568 E-01 , 

3 

y = I b.xj , 500°E < T < 1500°E (B.12h) 

i=o 

where h^ = 0.18501782 E 00, = -0.59566059 E-01, 

bg = 0.74697853 E-01, b^ = -0.20197244 B-01 , 

3 

y = I , 1500°E < T < 3000°E (B.12c) 

i=o 

where = 0.1 1932987 E 00, °1 = 0.74511450 E-01 ,■ 

Cg = -0.18387187 E-01, = 0.17363686 E-02, 

3 

y = I d.xj , 3000°E < T <. 6000 (B.12d) 

i=o 

where = 0.16704180 E 00, d^ = 0.29563507 E-01, 

d2 = -0.41591628 E-02, d^ = -0.22257059 E-03. 

The number of data points used in deriving the above 
relations, the maximum error involved at the data points, ard 
the sum of squared errors at the data points in the various 
cases are given in Table B.1. Though these relations are 
obtained for foot -pound -second system of units, these can be 
used for M.K.S. and S.I. units by using suitable multiplication 


factors . 
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APPENDIX C 

EXPRESSIONS FOE [AK], [BKj, . . [ DM J 


The following notation is used for convenience; 


u. 

1 

0 

dz 1 i = 1 , 2 , 

• • . , n 

(C.1) 

^i 

— 

1 

*H 

H 

il 

? i = 1 , 2 , 

• . . , n 

(0.2) 

^i 

= 

0 

cos^[ (©2 - ) f 

+ ] dz ; 




i = 1, 2, 

. . . , n 

(0.3) 

s, 

1 

1 

r 

= J z 

0 

^ sin 2 [ (©2 - 6-| ) 

Y + ] dz 1 




i = 1, 2, 

. . • 9 n 

(0.4) 

where 

and ©2 denote the values of 

pre twist at nodes 1 

and 2 

respectively of the 

element . 




As the nature of w^, , v^ and v^ is same except for 

their positions in the stiffness and mass matrices, one can use 
w to denote any one of the quantities w^, , v^ or v^ and in a 

similar manner the set (u^ , Ug s can be used to repre- 

sent anyone of the sets (u^ , U 2 » u^ , u^), (u^ , Ug , Urj, Ug) , 

(ug j ^10’ ^11 ’ ^12^ ^15’ ^16^* 
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w(z) = ^ (2z^ - 31z^ + 1^) - ^ (z5 _ 21z^ + l^z) 

+ (31z^ - 2z’) - ^ (z3 - 1^2) (C.5 

1 3- 

II = (6z^ - 61 z) - {3 z 2 _ + i2) + 1^1 

(61z - 6z^) - ^ (32^ - 2l2) (C.6) 

^ (I2z - 61) 1 (6z - 41) + -| (61 - 12z) 

dz 1^ i‘^ id 

- ^ (6z - 21) (C.7) 


By letting 1*1^^ (i = I,-.., 4; j = i,...,4? k = 1,...,7) 

denote the coefficient of z^~^ l"^”^ for u. u . term in the 

-2 ^ ^ 
expression of w , 

Q- • •]_ (4 “ 1^««9^49 3 “ dy»»»y45 k — 1j«»*y3) 

X j J 

k -I 5 -k - - 

the coefficient of z 1 for term in the expression of 

^dz^ ’ 


^ i_ (i ~ 1>***>45 D ~ d}«».j45 k — 1j».»j3) 

1 j j j ji 

k— 1 3— k - - 

the coefficient of z 1 for term in the expression of 

^2,1;: 2 
/dw\ 
v Q ' 1 
dz^ 
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= 1»*«*j 45 j the index coefficient of 1 to 

account for the difference in index of 1 due to multiplication 
of rotational degrees of freedom u^ and u^ the displacement 
degrees of freedom u, and u. , the values of P. , Q. . , , 

^ and H. . can he obtained as shown in Tables C.1 and C.2. 

-L J J J- , j 


Evaluation of [rk] : 

As the procedure for the derivation of [4K], [BK J, .. 
[dm] is same the expression f or [ BKJ is derived here as an 
illustration. 

T 1 3 v-i 2 

[Uq u^^ u^ 2] [BK][Ug u^Q u^^ ^^12]= / 

0 3 Z 


1 .2- 2 

J El a- 

o az 


(C.8) 


where w = 


v^ and 


u. 


Ur 


Ur 


u, 


Ur 






u 


10 


u 


u 


11 


12 


Putting the value of I and w in the Eq. (C,8), 

«y 



218 


1 ^2- 2 
; Ei™{^) dz 
0 az 


/ + (ly.y. - 


OOB^ ( {®2 " ®1 ^ f t ^ (ISz - 61) - ^ (6z - 41) 


+ — (61 - 12z) - -| (6z - 21)] 2 dz 


(C.9) 


with 


'1,1 ~ Coefficient of = Coefficient of Ug Ug 

H 1 

= ^To ^ / [ <(a.z^ + a„lz^ + a^l^z^ + + ap-1'^) 

121 ^ 1 2 3 4 5 

+ { (d^z‘^ + d2lz^ + d^l^z^ + d^l^z + d^l"^) 

- (a^z^ + a2lz^ + a^l^z^ + a^l^z + a^l^)} 


COS { (©2 - ^ 


^^1,1,1 ^ ®1 ,1 ,2 "*'^1,1,3 


(Hi^1+ i) ^i^®1,1,1 ^8-i ^ ^i+1^1,1,2^7-i 


^i+2 ^1,1,3^6-i^ ^*^i “ ^i^ ^ ^i^l ,1 ,1^8-j 

^i+1^1,1 ,2^7-i ^ ^i+2^1,1,3^6-i^ 
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5 3 






(di aj_) 


(C.10) 


This relation can be generalized as 


5 3 


121’° ® ^ 


1) ,1, I 


+ (dj^ - a^) «I(i+j.i)T^(9_i.j)El,j,3>! 


I zr J ~ 


10 ^ 

101 i-1 D 




U „ . . R ^ T ^ > + (‘^-5 ” ^ 

Q-x-n I.J.l 1 1 


"V/n ,• J ^I^T T 5 I = 1 , . . . ,4; J - I ? • • • >4 


(i+j+Hj^j) (9-i-D) I>J»D 


(C.ll) 


Similarly 


5 3 




+ (a^ - d.) {E(i+3+H ,)’'9-l-j®I,J,3 ’ 


1 = 1 ^ • ^4 j 1 — Ij*** j4 


(G.12) 
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3 5 




I = *ly*«*j4‘5 ^ 


(0.13) 


E 


5 3 


121 


10 


.L .L ^(i+D+% j)^(9-i-3)^I,J»d^’ 

1=1 :]=i ^9^ 


I = 1 j . ♦ • fA 5 J — !)•#• >4 


(C.14) 


3 T 

-f .1 I f^°i^(i,j+H3. )^(l1-i-3)^I,J,3^ ’ 

1 1=1 .1=1 *^9^ 


I =r 1,.*.^45 J — I,*..j4 


(0.15) 


"m 


121 


10 


5 5 

i=1 .1=1 


+ ^ ^(i+3+Hj^j)'^(11-i-D)'^I»J»3^ ^ 


I = 1,...,4| J = I,..*54- 


( 0 . 16 ) 


Pm 5 5 


ill ll [ j ' 


121*'° i=1 3 


- aj^) <I'(i+3+Hj^j)^(11-i-j)®I,J,3^^’ 


I = 1,...,45 j - 1»***}4- 


(0.17) 
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Prr. 5 5 

JD V V 




I — 1,...}4| J — 


(C.18) 


Evaluat ion of [ EK] and [ EK] ; 


1 3w, 2 

‘ -n/" _ N f D \ 


[u^ U 2 u^] [EK] [u^ U 2 u^] = / P(z)(-^) dz 


/ -P^A a2 [ (eL + - ez^ - z^) 


(e + Zg)z - -gz^] 


3W, 2 

(-^) dz 


/ p A £^2 (eL +^l 2 - ez^ - 4 ) dz 


/ p^A (e + Zg)z (— ) dz - / I £ 12^2 


(C.19) 


Thus 


p_ £1 12 


iL ii ■ 


°(9-i-j)®I,J,3 ^ 


3 5 


.T . -i 1 


(e + Zg) [°i^(i+3+Hj^j)^(10-i-3) I,J»D 
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P 

m 



ill jl-i 




(C.20) 


Similarly 


2 P Q 
m 


3 7 


i ^I’■T T 4 3 ? 




I = 1 , ... ,4; J = I, ... >4 


(C.21) 
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TABLE C.1 



1 

1 

0 

144.0 

- 144.0 

36.0 

36.0 

- 72.0 

36.0 

0.0 

0.0 

1 

2 

0 

- 144.0 

144 .0 

- 36.0 

- 36.0 

72.0 

- 36.0 

0.0 

0.0 

1 

3 

1 

- 72.0 

84.0 

- 24.0 

- 18.0 

42.0 

- 30.0 

6.0 

0.0 

1 

4 

1 

- 72.0 

60.0 

- 12.0 

- 18.0 

30.0 

- 12.0 

0.0 

0.0 

2 

2 

0 

144.0 

- 144.0 

36.0 

36.0 

- 72.0 

36.0 

0.0 

0.0 

2 

3 

1 

72.0 

-84 .0 

24.0 

18.0 

- 42.0 

30.0 

- 6.0 

0.0 

2 

4 

1 

72.0 

- 60.0 

12.0 

18.0 

- 30.0 

12.0 

0.0 

0.0 

3 

3 

2 

36.0 

- 48.0 

16.0 

9.0 

- 24.0 

22.0 

- 8.0 

1,0 

3 

4 

2 

36.0 

- 36.0 

8.0 

9.0 

- 18.0 

11.0 

- 2.0 

0.0 

4 

4 

2 

36.0 

- 24.0 

4.0 

9.0 

- 12.0 

4.0 

0.0 

0.0 
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Elements 


where 

a-| 

^12 

^13 ^ 

^14 ^ 

®-22 

^23 "" 

^24 ^ 

^•33 " 


APPENDIX B 


ELEMENTS OP MATRICES [A], [b], [C] AND [ D ] 


of the Matrix [ A ] s 


[A] 


1200 


ail 

ai2 

^13 

3-14 


^22 

^23 

^24 



^33 

s-34 


Symmetric a^^ 


(D.1 ) 


3 (396 a^ + 441 a2 + 504 + 630 a^ + 1260 a^) 

■ a^i 

. (114 a^ + 126 ag + 147 a3 + 210 a^ + 630 a^) 

. (282 a^ + 315 a2 + 357 a3 + 420 a^ + 630 a^) 

'11 

■ ^"13 

■ ^14 

- (36 a^ + 42 a2 + 56 a3 + 105 a^ + 420 a^) 
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®-34 ^ 1 + 84 ag + 91 + 105 a^ + 210 a^) 

= Y (204 a^ + 231 ag + 266 a^ + 315 a^ + 420 a^) 

Elements of Matrix [b] ; 



N 1 

= 

i (432 

+ 756 + 1512 

bi2 

= 

- b^^ 


bi3 

= 

- (90 

+ 126 Cg + 126 C 3 ) 

^14 

= 

- (36 

+ 0 + 126 c_) 

3 

^22 

= 

^11 


^23 

= 

-^3 


^24 

= 

-"14 


b33 

=: 

1(24 + 

42 C 2 + 168 c^) 

^34 

= 

- 1(18 

+21 + 42 c^) 

\4 


1(108 

+ 126 + 168 c^) 
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Elements of Matrix [C ] ; 






C11 c^2 

"13 

"14 



1 — 1 

0 

11 

o^ B 

0 

°22 

°23 

C33 

°24 

"34 





Symmetric 


C44 


= 

1(38 

+ 108 

c^ + 468 c^) 



^2 

= 

1(46 

+ 81 c 

2 + 162 Cj ) 



°13 

= 


+ 21 Cg + 66 Qj) 



°14 

= 

( 13 . c 

1 i 2 

+ 19 

=2 + 39 03) 



C\J 

C\J 

o 

= 

1(290 

+ 360 

C2 + 468 C^) 



Q 

ro 

=: 

- l2(^ 
J- V 2 

Cl +.21 Cg + 39 c^) 



CM 

O 


i^r^ c 

-L l 2 1 

+ 45 

C2 + 66 c^) 



C33 

= 

1^(2 0, 

+ -2 c 
+ 2 °2 

+ 12 c^) 



°34 

= 

-l5(f c 

+ ^ 
12 

c 2 + 9 ^3 ^ 



*^44 

= 

1^(5 0, 

+ ^c 

2+12 c^) 




(D.3) 
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Elements of Matrix [d]; 







dll ^^12 *^13 ^14 



p 

^22 *^23 ^^24 


[d] = 

m 



1260 

d33 d^^ 




<^44 




u_ U 

d^1 

= i (15 a 

+ i5 . 
•1 + 2 ' 

3-2 + 36 a^ + 63 a^ + 126 a. 

di2 

= -111 



^13 

- -(-1^ a 

^41 

u. 21 

+ — a. 

+ Jl,+21a +ma) 

2^^ 2 ^4^ 2 0 

^14 


+ -^ a 

4 ®'2 

+ 3 + 0 - ^ a^) 

*^22 

= d^^ 



d23 

= - d^^ 



C\J 

= -^14 



'^33 

= l(a^ + 

1 1 

^ ^2 + 

2 .a3 + 2 a^ -•- H a^) 

S4 

= -i(|« 

1 ^ ‘ 

^2 + I ■" 4 ^4 2 

<^44 

= 1 ("^ 


P + 9 a^ + a^ + 14 a^) 


(D.4) 
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APPEiroiX E 

EESCRIPTIOE OP THE COMPUTER PROGRiiMME 


The computer programme for the automated optimum design 
of axial flow gas turbine stage is written in EORTRilJ IV 
language and it consists of 21 subroutines. Most of the 
computations are made on IBM 7044 at Indian Institute of 
Technology, Kanpur, but some of the computations are also done 
on IBM 360 of Delhi University. The purpose of the various 
subroutines is given below s 

(1) ALOAD: To formulate the load vectors due to pressure and 
gas bending toxoes in the directions of the degrees of 
freedom of a turbine blade . 

(2) CBMD: To perform the Choleshy decomposition of symmetric 
banded matrices. This routine stores only the upper 
triangular band of the decomposed matrix and the diagonal 
of the matrix is stored in the first column. 

(3) ZOLP: To solve the equilibrium equations from the upper 
triangular band of the decomposed stiffness matrix. 

(4) MATINV; To find the inverse of a real square matrix. 

(5) EIGEU: To compute the eigen values and eigen vectors of 
the generalized Ritz problem by power method. 
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(6) BASICS To obtain the partial eigen solution of a structure 
by Rayleigh -Hit z sub-space iteration technique . It calls 
ZOLP, MATIirV and EIGEN for solving the generalized Ritz 
problem . 

(7) STRESS s To find stresses, deflections and eigen values 
of rotating, tapered and twisted beams of rectangular 
cross-section. It calls the subroutines ABO AD, CBANI), 

ZOIE and BASIC. 

(8) EIUs To perform the Cholesky decomposition of symmetric 
matrices. This will not take advantage of the banded 
nature of the matrices. 

(g) SOlVEBi To solve a set of simultaneous equations using 

the lower and upper triangular matrices obtained from the 
subroutine ELU. This will not take advantage of the banded 
nature of the original coefficient matrix. 

(10) CAMJl; To obtain the coefficients of the square and cube 
of a fourth order polynomial. 

(11 ) CAROL: To evaluate the value of a n^ order polynomial 
when the values of the coefficients and the variable are 

given . 

(12) POLEIT: To fit a polynomial for a set of data points 
using least squares method. It calls ELU and S0L7EB 

routines . 
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(13) VINT: To evaluate the integral of a general n"^^ order 
polynomial in terms of the coefficients of the polynomial 
and the limits of integration. 

( 14 ) SECT: To obtain an eq,uivalent rectargular cross-section 
of an airfoil section. 

( 15 ) YPLOSS: To calculate the profile loss coefficients for 
the given blading . 

(16) OPT: To obtain the value of specific heat at constant 
pressure for a given temperature range. 

( 17 ) ETN: To evaluate the objective function penalty function 
and the constraints of the optimization problem for any 
given design vector. It calls STEESS, SECT, YPLOSS and 
OPT subroutines. 

(18) GRAIN: To calculate the gradients of objective function, 
penalty function and constraints using backward difference 
method. It calls the subroutine ETN as many times as there 
are design variables . 

( 19 ) S10PE2; To evaluate the slope of the penalty function. 

It calls GRALN subroutine . 

(20) STEPX3; To implement the one dimensional search method by 
cubic interpolation. It calls the subroutines ETN and 


SLOPE 2 . 
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(21) s To implement the Davidon-Pletcher-Powell variahle 

metric method of unconstrained minimization. It calls PTIT, 
GRADM and STEPX3. 

The interrelation between the various subroutines is 
shown in the following chart. 


MIR 


MmM4(21 ) 


PTN(17) 


gr1dn(18) 

PTR(17) 


r 


PIIKIT) 


r 


I 


STEPX3{20) 

SL0PE2(19) 


GRAPRde) 


, , 

YPL0SS(15) SECT(14) CPT(16) STEESS(7) 

r \ ^ ^ ^ 

ELU(8) S0LVEB(9) POLEIT(12) CMUIi(lO) VIRT(13) 


ELU(8) SOLVES (9) 


C'APOL(II) 


1 \ I . 

ALOAD(I) CBMD(2) Z0LP(3) BApC(6) 


Z0LE(3) 


"T 

MATI]W(4) 




EIGER(5) 
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[K] 

[M] 

[-X] 

P 

[Q] 


= generalised stiffness matrix 

= generalised mass matrix 

= initial approximate vectors 

= degree of freedom 

= number of required first few modes 

= maximum number of subspace iterations allowed 

= eigen vectors of the projjected operators 


Enter with [k]J!m],[x] jU^jP 


Set ITER = 0 


' 9 0 ) 


Set ITER = ITER +1 

iteration counter for sub-space iteration. 


IS ITER > TIER. 


max 



Yes 



I Solve for [Y] 
=[M] X 


Compute [K] and [M] 

[s] = pf [^[■y] 

[M] = [Y]^ [m3[y] 


Divide jth row of [g], [M] and jth column 
of for j = 1, P 


continued. . . 
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Pig. B.1 Rayleigh-Ritz Sub -space Iteration Algorithm for 
Determining Eigen Solutions 






